Discrete port-Hamiltonian system model of a single-reed woodwind instrument

Time-domain simulation of woodwind instruments typically involves the development of separate discrete-time sub-models for the excitation mechanism and the resonator. These components have largely been modeled via digital waveguide or finite-difference time-domain (FDTD) methods. We present a separate approach based on the modular and energy-based port-Hamiltonian system (PHS) framework. We recast the three main components of a woodwind instrument—the single-reed, the bore, and the tonehole—as PHS models and incorporate novel elements in each derivation. In the beating reed model, we make use of recent work on energy quadratization to formulate a linearly implicit scheme of the nonlinear Hunt-Crossley contact force coupled to a nonlinear Bernoulli flow. In the horn model, we discretize a distributed PHS representing the horn equation with a generalized symplectic Störmer-Verlet scheme, verifying previously proposed FDTD schemes. In the tonehole model, we propose a new low-frequency model of the tonehole and model note transitions with a switching PHS. The benefit of describing each element as a PHS is demonstrated by the ability to interconnect all sub-models in a modular and energy-conserving manner to simulate a complete instrument. Simulations are performed on a test instrument and the numerical stability of the overall scheme is demonstrated.

1 Introduction

Research into time-domain woodwind physical modeling synthesis began to coalesce with the publication of McIntyre, Schumacher, and Woodhouse’s seminal 1983 article (McIntyre et al., 1983). In their article, the authors demonstrate an efficient time-domain approach for modeling different musical instruments including the clarinet. The article also presents concepts that are found in physical modeling synthesis research to this day, namely, the subdivision of an instrument into modular components, including a nonlinear excitation model and a passive resonator model. In a woodwind instrument, these components correspond to the reed and the instrument bore, respectively. The methods proposed in McIntyre et al. (1983) were further developed in the field of computer music by digital waveguide (DWG) synthesis (Smith, 1986; Smith, 1992).

DWG synthesis is one of the most prevalent methods for modeling woodwind and brass instruments (Välimäki, 1995; Scavone, 1997; van Walstijn and Campbell, 2003; Mignot et al., 2010). DWG synthesis is based on a discretization of the traveling wave solution of the one-dimensional wave equation. The result is a digital system composed of delay lines—for propagating the traveling waves—and filters, for reflecting waves, modeling losses, and approximating fractional delay lengths. Reed dynamics are incorporated by modeling the reed-tip as a mass-spring-damper system (Scavone and Cook, 1998). The popularity of the DWG method can be attributed to its efficiency and the fact that the DWG models consist of delay lines and filters that fit naturally into the framework of digital signal processing (DSP). Another benefit of the DWG method is its modularity. A DWG model can be viewed as an interconnection of simpler component sub-models.

More recently, energy-based numerical simulation techniques using finite-difference time-domain (FDTD) methods have become a popular avenue for physical modeling synthesis (Bilbao, 2009). FDTD methods directly discretize the partial differential equations (PDEs) governing the physical model. In contrast, DWG methods rely on a discretization of analytical solutions to a system—as in the case of modeling wave propagation. This, however, typically limits the DWG method to modeling cylindrical or conical bore profiles with further limitations on the stability of the method for convex bore shapes (Berners, 1999). An approach to modeling bores with profiles that are $C^2$ smooth is presented in Mignot et al. (2008). FDTD methods have been used to model the single-reed (Bilbao, 2008), bores with varying cross-sections (Bilbao and Chick, 2013), and viscothermal losses (Bilbao et al., 2015a; Bilbao and Harrison, 2016). FDTD methods have been used to model brass instruments as in Harrison-Harsley (2018); Willemsen (2021). Although FDTD methods are more computationally expensive than scattering methods for simple one-dimensional systems, they are more general and flexible making it possible to simulate non-trivial bore profiles which would not be stable in a DWG implementation. Given the inherent and desirable nonlinear behavior of musical instruments, energy-based techniques such as FDTD methods provide a general framework for the analysis of the stability of nonlinear numerical schemes. In comparison, analysis of scattering-based methods is limited by linear time-invariant (LTI) systems theory.

An alternative general energy-based modeling approach—and the focus of this article—is the port-Hamiltonian system (PHS) framework introduced in Maschke and van der Schaft (1993). PHSs combine energy-based analysis techniques within the Lagrangian and Hamiltonian framework from classical mechanics with the network modeling framework from electrical engineering (van der Schaft, 2006). The PHS framework provides a helpful procedure for the design and analysis of a continuous-time physical model as energy conservation is baked into the system formulation and models can be composed in a modular fashion (Duindam et al., 2009). PHSs have been primarily used to model nonlinear audio circuits in the related field of virtual analog (Falaize and Hélie, 2016a). Regarding physical modeling synthesis, PHSs have been used to model nonlinear strings (Hélie and Roze, 2016), electric pianos (Falaize and Hélie, 2017b), the lip-reed mechanism in brass instruments (Lopes and Hélie, 2016; Lopes, 2016), and a mobile vocal tract (Wetzel et al., 2019). While the PHS framework seems ideal for sound synthesis—the method is modular, energy-based, and formulated to be energy conserving—it is not without its drawbacks. To guarantee energy conservation in discrete PHSs, the primary discretization method is the discrete gradient method (Yalçin et al., 2015; Falaize and Hélie, 2016a). Higher-order methods for discretizing a PHS are discussed in Müller (2021). The discrete gradient method results in numerical dispersion, i.e., frequency warping, akin to the bilinear transform (Harrison-Harsley, 2018). For musical applications, this is undesirable as frequency warping can affect the tuning of simulated instruments. It is possible to minimize numerical dispersion through the careful design of explicit numerical schemes which are instead conditionally stable. Additionally, discrete PHSs in physical modelling synthesis are often relegated to lumped system modeling. Discrete distributed systems for modeling wave propagation have been handled using DWG (Lopes and Hélie, 2016) or modal decomposition techniques (Hélie and Roze, 2016; Falaize and Hélie, 2017b). However, modal decomposition is only valid for certain boundary conditions and if a decomposition exists.

The primary goal of this article is to present a PHS description of the three main components of a woodwind instrument: a lumped single-reed excitation mechanism, a one-dimensional bore model with variable cross-section based on the horn equation, and a lumped model of a tonehole. The secondary goal of this article is to draw similarities between the PHS and FDTD literature, presenting PHSs in a way that will be familiar to practitioners in the field of physical modeling synthesis and discretizing them with FDTD methods. Similar models of the reed, bore, and tonehole have been proposed based on FDTD methods in Chatziioannou et al. (2019); Bilbao and Harrison (2016) and DWG methods in van Walstijn and Scavone (2000), respectively and this article contributes refinements to these models. We include recent research on energy quadratization methods (Ducceschi et al., 2021; van Walstijn et al., 2024b) to design a linearly implicit scheme for the Hunt-Crossley contact in the reed system and propose a new low-frequency model of the tonehole that better approximates frequency-domain models in the literature. Energy conservation during note transistions is handled by modeling the tonehole as a switching PHS. We describe the lossy horn equation model proposed by Bilbao and Harrison (2016) as a distributed PHS and discretize the lossless wave propagation using the structure-preserving symplectic Störmer-Verlet method (Hairer et al., 2000), affirming the scheme used in their article. Furthermore, by characterizing each individual sub-model as a PHS we are able to interconnect each element in an energy-conserving and modular fashion and develop a system for modeling bores with arbitrary geometries and tonehole placements.

This paper is organized as follows: Section 2 will introduce port-Hamiltonian systems. Section 3 will discuss FDTD discretization methods and introduce the discrete gradient and Störmer-Verlet methods as FDTD methods. Section 4 provides an example of modeling and discretizing an RLC system model in the PHS framework. Section 5 will discuss the single-reed excitation including contact dynamics with the mouthpiece. Section 6 will present a distributed PHS model of the woodwind bore and Section 7 considers a model of the woodwind tonehole based on a switching PHS formulation. Section 8 will display simulation results with Section 9 concluding the article.

2 Port-Hamiltonian systems

Port-Hamiltonian systems (PHSs), introduced in Maschke and van der Schaft (1993), approach physical modeling through the Hamiltonian equations of motion and combine these equations with network theory. The formalism provides a framework for a geometric description of models based on the interconnection of sub-systems through power-conserving interconnections (van der Schaft, 2006). As a consequence, energy conservation is baked into the PHS formulation through preservation of a power balance. Energy-conserving numerical methods—if discretized properly—will preserve passivity and avoid numerical instability. The PHS approach uses energy as a lingua franca between domains, making the framework well-suited for multi-physical simulations. Musical instruments are inherently multi-physical as the mechanical actions of the player are transformed into acoustic energy. Musical instruments are also inherently nonlinear. Nonlinearities such as the reed contact with the mouthpiece or the slip-stick action of a bowed string directly affect the perceived timbre of an instrument during performance (Chaigne and Kergomard, 2016). These reasons naturally motivate the use of PHSs for modelling musical systems. In the context of physical modeling, the PHS formalism is attractive because it provides a general form that can be used to characterize a variety of physical systems. This continuous-time formulation can then be discretized in a structured manner to produce a discrete-time numerical scheme.

In this section, we will review Hamiltonian dynamics and then the form of finite-dimensional PHSs. Theoretical specifics such as the definition of a Dirac structure are left out of this review. Instead, we circularly define Dirac structures as structures that preserve the power balance and are a consequence of Kirchhoff’s current and voltage laws. Interested readers are referred to the texts (van der Schaft, 2006; Duindam et al., 2009) and the tutorial (Hélie, 2022).

2.1 Hamiltonian dynamics

The Hamiltonian equations of motion for a system with $N$ particles consist of $2N$ first-order equations

${\partial }_t{q}_n={\partial }_{p_n}\mathcal{H},(\mathrm{1a})$

${\partial }_t{p}_n=-{\partial }_{q_n}\mathcal{H},(\mathrm{1b})$

where ${\partial }_x$ is the partial derivative with respect to a general variable $x$. Each particle in the system has positions $\left[q_1,q_2,\dots ,q_N\right]$ and momenta $\left[p_1,p_2,\dots ,p_n\right]$ (Goldstein, 1980). $\mathcal{H}$ defines the Hamiltonian—the total energy of the system

$\mathcal{H}\left(\mathbf{p},\mathit{q}\right)={\displaystyle \sum _n}\mathcal{T}\left(p_n\right)+\mathcal{V}\left(q_n\right),(2)$

where $\mathcal{T}$ is the component defining kinetic energy and $\mathcal{V}$ is the component defining potential energy. These components are dependent on the momentum and position of each particle. The system described by the Hamiltonian equations is not general. Notably, damping is omitted. This motivates the use of the more general PHS framework.

2.2 Port-Hamiltonian systems

Port-Hamiltonian systems describe the interconnection of energy storage elements, energy dissipating elements, and power conserving elements using general power-conjugate variables of effort, $\mathbf{e}$, and flow, $\mathbf{f}$ (Duindam et al., 2009). The product of effort and flow variables is equal to the power of the system, that is, the change in the total energy over time $({\partial }_{t}H)$. Effort and flow can be given physical meaning in electrical systems as voltage and current (Section 4), in mechanical systems as force and velocity (Section 5), and in acoustic systems as pressure and volume velocity (Section 6). The effort and flow variables are related to one another by a set of internal state variables which is the integral of the effort with respect to time. Taking the system in Equations 1a, 1b as an example, the state variables are the positions and momenta of each particle. The efforts are the gradient of the Hamiltonian with respect to the state variables and the flow is the time-derivative of each state variable. By the chain rule, the product of the effort and flows corresponds to power.

A useful representation of a dynamical PHS with state $\mathbf{x}$, input $\mathbf{u}$, and output $\mathbf{y}$ is the input-state-output port-Hamiltonian system with direct feed-through (Duindam et al., 2009)

$\underset{\mathbf{f}\left(\mathbf{x}\right)}{\underbrace{{\partial }_t\mathbf{x}}}=\underset{\mathbf{A}}{\underbrace{\left(\mathbf{J}-\mathbf{R}\right)}}\underset{\mathbf{e}\left(\mathbf{x}\right)}{\underbrace{\nabla \mathcal{H}\left(\mathbf{x}\right)}}+\underset{\mathbf{B}}{\underbrace{\left(\mathbf{G}-\mathbf{P}\right)}}\mathbf{u},(\mathrm{3a})$

$\mathbf{y}=\underset{\mathbf{C}}{\underbrace{{\left(\mathbf{G}+\mathbf{P}\right)}^T}}\nabla \mathcal{H}\left(\mathbf{x}\right)+\underset{\mathbf{D}}{\underbrace{\left(\mathbf{M}+\mathbf{S}\right)}}\mathbf{u}.(\mathrm{3b})$

$\mathbf{e}$ and $\mathbf{f}$ are the state dependent efforts and flows. $\nabla$ is the gradient operator and $\mathbf{J}$ is a skew-symmetric matrix ($\mathbf{J}$ = $-{\mathbf{J}}^T$) representing the energy-preserving interconnection between storage components. $\mathbf{R}$ is a matrix representing dissipative components and must be positive-definite in a passive system. $\mathbf{G}$ encapsulates the relationship between external ports and the internal states with $\mathbf{P}$ representing the dissipation between external ports and the storage components. $\mathbf{M}$, also a skew-symmetric matrix, is the energy-preserving interconnection and $\mathbf{S}$ the dissipation between the input and output ports. The resulting PHS structure consists of a mixed set of differential and algebraic equations (DAE) (Duindam et al., 2009). The system in Equations 3a, 3b will also be familiar as the state-space equations from dynamical systems theory and signal processing (Kailath, 1980; Smith, 2010). Equation 3a is the state update equation and Equation 3b is the output equation. If damping, input, and output variables are omitted, and we define $\mathbf{x}={\left[\mathit{q}\mathbf{p}\right]}^T$, $\mathbf{f}(\mathbf{x})={\partial }_t\mathbf{x}$, and $\mathbf{e}(\mathbf{x})=\nabla H(\mathbf{x})$ then Equation 3a is equivalent to Equation 1b.

Defining the total effort $\mathbf{E}=\left[\nabla \mathcal{H}\left(\mathbf{x}\right)\mathbf{u}\right]$, and flow variables $\mathbf{F}=\left[{\partial }_t\mathbf{x}-\mathbf{y}\right]$, we can rewrite Equation 3 as

$\mathbf{F}=\underset{\mathcal{J}}{\underbrace{\left[\begin{array}{cc}\hfill \mathbf{J}\hfill & \hfill \mathbf{G}\hfill \\ \hfill -{\mathbf{G}}^T\hfill & \hfill -\mathbf{M}\hfill \end{array}\right]}}\mathbf{E}-\underset{\mathcal{R}}{\underbrace{\left[\begin{array}{cc}\hfill \mathbf{R}\hfill & \hfill \mathbf{P}\hfill \\ \hfill {\mathbf{P}}^T\hfill & \hfill \mathbf{S}\hfill \end{array}\right]}}\mathbf{E},(4)$

where we have the block skew-symmetric matrix $\mathcal{J}$ and the block dissipative matrix $\mathcal{R}$. Due to the property for any skew-symmetric matrix ${\mathbf{x}}^T\mathbf{J}\mathbf{x}=0$, the inner product of efforts and flows expresses the power balance

$⟨\mathbf{E},\mathbf{F}⟩=\underset{\begin{array}{c}={\partial }_t\mathcal{H}\\ \text{Energy\hspace{0.17em}variation}\end{array}}{\underbrace{\nabla {\mathcal{H}}^T\left(\mathbf{x}\right){\partial }_t\mathbf{x}}}=\underset{\begin{array}{c}\mathcal{B}\\ \text{External\hspace{0.17em}power}\end{array}}{\underbrace{{\mathbf{u}}^T\mathbf{y}}}-\underset{\begin{array}{c}\mathcal{Q}>0\\ \text{Dissipated\hspace{0.17em}power}\end{array}}{\underbrace{{\mathbf{E}}^T\left[\begin{array}{cc}\hfill \mathbf{R}\hfill & \hfill \mathbf{P}\hfill \\ \hfill {\mathbf{P}}^T\hfill & \hfill \mathbf{S}\hfill \end{array}\right]\mathbf{E}}}.(5)$

Equation 5 states that the variation in the total stored energy is equal to the energy flowing in and out of the system, $\mathcal{B}$, minus the energy that is dissipated over the dissipative elements, $\mathcal{Q}$. The system is passive if the change in stored energy is less than the energy flowing into and out of the system ${\partial }_{t}H\le {\mathbf{y}}^T\mathbf{u}$. This holds true if the dissipative elements satisfy

$\mathcal{R}=\left[\begin{array}{cc}\hfill \mathbf{R}\hfill & \hfill \mathbf{P}\hfill \\ \hfill {\mathbf{P}}^T\hfill & \hfill \mathbf{S}\hfill \end{array}\right]\ge 0.(6)$

We will also define a conjugate PHS where the variables ${\partial }_t\mathbf{x}$ and $\nabla \mathcal{H}\left(\mathbf{x}\right)$ are interchanged. The total effort will then be $\mathbf{E}=\left[{\partial }_t\mathbf{x}\mathbf{u}\right]$ and the total flow will be $\mathbf{F}=\left[\nabla \mathcal{H}\left(\mathbf{x}\right)-\mathbf{y}\right]$. The conjugate can be seen as a particular case of a hybrid PHS (van der Schaft and Jeltsema, 2014, Chapter 5) where all the storage efforts and flows have been interchanged. This form will be seen in Section 6.2. A hybrid PHS, where one set of efforts and flows have been interchanged, is used to model the switching PHS tonehole in Section 7. These alternative representations still maintain the power balance in Equation 5 as each row in the total effort and flow corresponding to storage elements is still a set of power-conjugated variables. Hybrid Dirac structures are used in Müller and Hélie (2021) along with dedicated solvers based on projection methods and a Newton iteration.

An important class of Hamiltonian functions are quadratic Hamiltonians. Quadratic Hamiltonian functions appear frequently in the systems we are interested in and have the following form and gradient,

$\mathcal{H}\left(\mathbf{x}\right)=\frac{1}{2}{\mathbf{x}}^T\mathbf{L}\mathbf{x},\hspace{1em}\nabla \mathcal{H}\left(\mathbf{x}\right)=\mathbf{L}\mathbf{x},(7)$

where $\mathbf{L}$ is a symmetric matrix. Passivity is guaranteed if $\mathbf{L}\ge 0$. Even if the Hamiltonian is not quadratic it can be turned into a quadratic function by a change of variables using energy quadratization methods (Lopes et al., 2015; Yang, 2016).

3 Discretization

Discretization methods are at the heart of any physical modeling synthesis algorithm and their design is of great importance. We will begin this section by first reviewing finite-difference time-domain (FDTD) methods, adopting the notation used in Bilbao (2009) for the rest of this article. We will first present the elementary discretization operators used in the design of numerical schemes and some useful identities given by Bilbao (2009). Then, we will present the two discretization methods used in this article. The first is the discrete gradient method, commonly used in PHS modeling, and the second is the symplectic Störmer-Verlet method. These methods are examples of mechanical integrators that preserve either the discrete energy or the symplectic form of the continuous-time system (Wendlandt and Marsden, 1997). The symplectic form is equivalent to the orientated area in state space. Symplecticity is a characteristic property of Hamiltonian systems and symplectic integrators are used to maintain this property in a numerical scheme (Hairer et al., 2000; Sanz-Serna and Calvo, 1994). It has been observed from numerical experiments (Chatziioannou, 2019) and backward error analysis (Hairer et al., 2003) that symplectic schemes provide a good discretization of a system’s dynamics at the expense of the discrete energy not being an exact sampling of the continuous-time energy function. Instead, a modified discrete energy function is conserved. In the physical models proposed in this article, we will utilize the discrete gradient method in scenarios where there is nonlinear behavior, for easier control over the scheme’s stability, (Section 5), or when the behavior of the system results in minimal numerical dispersion (Sections 6.2, 7). We will utilize the Störmer-Verlet method where the system is passive, and numerical dispersion is a concern (Section 6).

3.1 Finite difference operators

We define the following general shift operators applied to a function of time, $t$, dimension $z$, and a general variable $s$

$e_{s+}u\left(t,z,s\right)=u\left(t,z,s+{\mathrm{\Delta }}_s\right)=u_l^n\left(s+{\mathrm{\Delta }}_s\right),e_{s-}u\left(t,z,s\right)=u\left(t,z,s-{\mathrm{\Delta }}_s\right)=u_l^n\left(s+{\mathrm{\Delta }}_s\right).(8)$

Because FDTD methods represent variables at regular discrete time, $t=n{\mathrm{\Delta }}_t$, and spatial coordinates, $z=l{\mathrm{\Delta }}_z$, a shorthand is used where the temporal and spatial coordinates are absorbed into a superscript and subscript, respectively. Here, we describe these operators on a general variable $s$ to highlight that the operators are not limited to just time and spatial dimensions.

Elemental finite difference operators, represented by the operator $\delta$ approximate continuous-time partial derivatives ${\partial }_s$. Some elementary operators are the forward $({\delta }_{s+})$, backward $({\delta }_{s-})$, and centered $({\delta }_{s\cdot })$ difference operators

${\partial }_s\approx \left\{\begin{array}{c}{\delta }_{s+}=\frac{1}{{\mathrm{\Delta }}_s}\left(e_{s+}-1\right)\hfill \\ {\delta }_{s-}=\frac{1}{{\mathrm{\Delta }}_s}\left(1-e_{s-}\right)\hfill \\ {\delta }_{s\cdot }=\frac{1}{2{\mathrm{\Delta }}_s}\left(e_{s+}-e_{s-}\right)\hfill \end{array}\right..(9)$

Of use are the forward $({\mu }_{s+})$, backward $({\mu }_{s-})$ and centered $({\mu }_{s\cdot })$, which approximate unity,

$1\approx \left\{\begin{array}{c}{\mu }_{s+}=\frac{1}{2}\left(e_{s+}+1\right)\hfill \\ {\mu }_{s-}=\frac{1}{2}\left(1+e_{s-}\right)\hfill \\ {\mu }_{s\cdot }=\frac{1}{2}\left(e_{s+}+e_{s-}\right)\hfill \end{array}\right..(10)$

These operators are useful for collocating variables in discrete time. Higher-order derivatives can be approximated by combining these elementary operators.

3.1.1 Helpful identities and inequalities

Working with a general function $u$, a useful identity used throughout this article is (Bilbao, 2009)

$\left({\delta }_{t\pm }u\right){\mu }_{t\pm }u={\delta }_{t\pm }\left(\frac{1}{2}{u}^2\right).(11)$

This product identity is the core of the discrete gradient method described in Section 3.2. We have the identities relating different operators

${\delta }_{t\cdot }u={\mu }_{t-}{\delta }_{t+}u={\mu }_{t+}{\delta }_{t-}u,(12)$

${\mu }_{t+}=1+\frac{{\mathrm{\Delta }}_t}{2}{\delta }_{t+},(13)$

${\mu }_{t\cdot }={\mathrm{\Delta }}_t{\delta }_{t\cdot }+e_{t-},(14)$

and the following product identities

$\left({\delta }_{t\cdot }u\right){\mu }_{t\cdot }u={\delta }_{t\cdot }\left(\frac{1}{2}{u}^2\right),(15)$

$\left({\delta }_{t\cdot }u\right)u={\delta }_{t+}\left(\frac{1}{2}u\left(e_{t-}u\right)\right),(16)$

and,

$u\left(e_{t-}u\right)={\left({\mu }_{t-}u\right)}^2-\frac{{\mathrm{\Delta }}_t^2}{4}{\left({\delta }_{t-}u\right)}^2.(17)$

The last identity produces the inequality:

$u\left(e_{t-}u\right)\ge -\frac{{\mathrm{\Delta }}_t^2}{4}\left({\delta }_{t-}u\right).(18)$

3.2 Discrete gradient method

The discrete gradient method is an energy-conserving discretization method for Hamiltonian systems. It involves the approximation of the gradient operator by a forward difference with respect to the state

$\nabla \mathcal{H}\left({\mathbf{x}}^n\right)\approx {\delta }_{\mathbf{x}+}\mathcal{H}\left(\mathbf{x}\right)=\left[\mathcal{H}\left({\mathbf{x}}^n+{\mathrm{\Delta }}_{\mathbf{x}}\right)-\mathcal{H}\left({\mathbf{x}}^n\right)\right]\oslash {\mathrm{\Delta }}_{\mathbf{x}},(19)$

where ${\delta }_{\mathbf{x}}+$ is the forward difference with respect to the state vector $\mathbf{x}$ and $\oslash$ defines the Hadamard or element-wise division between vectors. The key to the discrete gradient method is defining ${\mathrm{\Delta }}_{\mathbf{x}}$ as the difference between the state at successive time steps

${\mathrm{\Delta }}_{\mathbf{x}}={\mathbf{x}}^{n+1}-{\mathbf{x}}^n.(20)$

Defining the inner product as $⟨\mathbf{x},\mathbf{y}⟩={\mathbf{x}}^T\mathbf{y}$ for two vectors $\mathbf{x}$ and $\mathbf{y}\in {\mathbb{R}}^N$, the discrete variation in energy is equivalent to the inner product of the discrete gradient and the forward difference of the state in time

${\delta }_{t+}\mathfrak{h}=⟨{\delta }_{\mathbf{x}+}\mathcal{H}\left(\mathbf{x}\right),{\delta }_{t+}\mathbf{x}⟩,(21)$

and the discrete gradient method observes a discrete chain rule. The presented definition of the discrete gradient is only one particular definition and other definitions exist (Yalçin et al., 2015).

For a quadratic Hamiltonian, the discrete gradient approximates the continuous-time gradient in Equation 7 via the midpoint rule. This corresponds to the forward average operator

${\delta }_{\mathbf{x}+}\mathcal{H}\left({\mathbf{x}}^n\right)=\frac{1}{2}\mathbf{L}\left({\mathbf{x}}^{n+1}+{\mathbf{x}}^n\right)={\mu }_{t+}\mathbf{L}{\mathbf{x}}^n.(22)$

Based on Equation 11, the discrete energy is conserved as

$⟨{\delta }_{t+}{\mathbf{x}}^n,{\mu }_{t+},\mathbf{L}{\mathbf{x}}^n⟩={\delta }_{t+}\left(\frac{1}{2}{\left({\mathbf{x}}^n\right)}^T{\mathbf{L}\mathbf{x}}^n\right).(23)$

Any PHS discretized with the discrete gradient method produces an implicit and unconditionally stable scheme. However, if the Hamiltonian is quadratic and the resulting scheme is linearly implicit, an explicit update form can be derived via a linear system solution. To derive the scheme approximating Equations 3a, 3b, we approximate ${\partial }_t$ with the forward difference in time Equation 9 and aprapproximate $\nabla \mathcal{H}$ as described in Equation 22. The resulting linearly implicit scheme has the following explicit update and output equations

${\mathbf{x}}^{n+1}={\mathbf{A}}_d{\mathbf{x}}^n+{\mathbf{B}}_d{\mathbf{u}}^n,(\mathrm{24a})$

${\mathbf{y}}^n={\mathbf{C}}_d{\mathbf{x}}^n+{\mathbf{D}}_d{\mathbf{u}}^n,(\mathrm{24b})$

with

$\mathbf{H}={\left(\mathbf{I}-\frac{{\mathrm{\Delta }}_t}{2}\mathbf{A}\mathbf{L}\right)}^{-1},\hspace{1em}{\mathbf{A}}_d=\mathbf{H}\left(\mathbf{I}+\frac{{\mathrm{\Delta }}_t}{2}\mathbf{A}\mathbf{L}\right),\hspace{1em}{\mathbf{B}}_d={\mathrm{\Delta }}_t\cdot \mathbf{H}\mathbf{B},(\mathrm{25a})$

${\mathbf{C}}_d=\frac{1}{2}\mathbf{C}\mathbf{L}\left({\mathbf{A}}_d+\mathbf{I}\right),\hspace{1em}{\mathbf{D}}_d=\frac{1}{2}\mathbf{C}\mathbf{L}{\mathbf{B}}_d+\mathbf{D},(\mathrm{25b})$

defined based on the matrices in Equations 3a, 3b. Unless the system parameters are time-varying, the discrete state-space matrices (with subscript $d$) can be computed and stored prior to the simulation. In the case that we have a conjugate PHS system where ${\partial }_t\mathbf{x}$ and $\nabla \mathcal{H}$ are interchanged, the discrete gradient method results in the following scheme

$\mathbf{L}{\mu }_{t+}{\mathbf{x}}^n=\mathbf{A}{\delta }_{t+}{\mathbf{x}}^n+\mathbf{B}{\mathbf{u}}^n,(\mathrm{26a})$

$y^n=\mathbf{C}{\delta }_{t+}{\mathbf{x}}^n+\mathbf{D}{\mathbf{u}}^n,(\mathrm{26b})$

and the explicit update matrices are then

$\mathbf{H}={\left(\mathbf{A}-\frac{{\mathrm{\Delta }}_t}{2}\mathbf{L}\right)}^{-1},\hspace{1em}{\mathbf{A}}_d=\mathbf{H}\left(\mathbf{A}+\frac{{\mathrm{\Delta }}_t}{2}\mathbf{L}\right),\hspace{1em}{\mathbf{B}}_d=-{\mathrm{\Delta }}_t\cdot \mathbf{H}\mathbf{B},(\mathrm{27a})$

${\mathbf{C}}_d=\frac{1}{{\mathrm{\Delta }}_t}\mathbf{C}\left({\mathbf{A}}_d-\mathbf{I}\right),\hspace{1em}{\mathbf{D}}_d=\frac{1}{{\mathrm{\Delta }}_t}\mathbf{C}{\mathbf{B}}_d+\mathbf{D}.(\mathrm{27b})$

3.3 Störmer-Verlet method

Symplectic methods such as the Störmer-Verlet method are specifically used to discretize partitioned systems such as those describing Hamiltonian dynamics. Consider a general non-autonomous system defined by two variables ($q$, $p$) and the input variable $u$. $q$ and $p$ are both dependent on a variable $x$ and related to one another by functions $f$, $g$, and the partial derivative with respect to $x$

${\partial }_{x}p=f\left(q,p,u\right),{\partial }_{x}q=g\left(q,p,u\right).(28)$

Normally $p$ and $q$ correspond to momentum and position, and ${\partial }_{t}q=p/m$. This is not necessarily the case in the general system. Hairer et al. (2003) describe a Störmer-Verlet scheme for discretizing general autonomous partitioned systems. We propose an extension of the method to non-autonomous systems

${\delta }_{x-}{p}_{x+\frac{1}{2}{\mathrm{\Delta }}_x}=\frac{1}{2}\left[f\left(q_x,p_{x+\frac{1}{2}{\mathrm{\Delta }}_x},u_x\right)+f\left(q_x,p_{x-\frac{1}{2}{\mathrm{\Delta }}_x},u_x\right)\right],(\mathrm{29a})$

${\delta }_{x+}{q}_x=\frac{1}{2}\left[g\left(q_{x+{\mathrm{\Delta }}_x},p_{x+\frac{1}{2}{\mathrm{\Delta }}_x},u_{x+\frac{1}{2}{\mathrm{\Delta }}_x}\right)+g\left(q_x,p_{x+\frac{1}{2}{\mathrm{\Delta }}_x},u_{x+\frac{1}{2}{\mathrm{\Delta }}_x}\right)\right].(\mathrm{29b})$

For linearly damped autonomous systems, the method in Equations 29a, 29b results in a contraction of symplectic area. In Chatziioannou (2019) it was found for a mass-spring-damper system that the contraction of the Störmer-Verlet scheme is a Padé approximation to the continuous-time contraction. Determining if this property holds for the class of all linear non-autonomous PHSs discretized with the proposed method and if the proposed method is symplectic in the non-autonomous case is outside the scope of this article.

The Störmer-Verlet method involves placing one variable, here $p$, on an interleaved grid at indices $x+\frac{1}{2}{\mathrm{\Delta }}_x$ and keeping the other variable, here $q$, on the normal grid. Positions between the interleaved and normal grid are related by the forward or backward averaging operators. For instance, in the case of the input variable

$u_x={\mu }_{x-}{u}_{x+\frac{1}{2}{\mathrm{\Delta }}_x},u_{x+\frac{1}{2}{\mathrm{\Delta }}_x}={\mu }_{x+}{u}_x,(30)$

Equation 29a can be understood to be centered on the normal grid. As such, the right-hand-side of Equation 29a averages $f$ evaluated at $p_{x+\frac{1}{2}{\mathrm{\Delta }}_x}$ and $p_{x-\frac{1}{2}{\mathrm{\Delta }}_x}$ with $q$ and $u$ fixed at $x$. Conversely, Equation 29b is centered on the interleaved grid. The traditional three-step form of the Störmer-Verlet scheme is given for the general method as

$p_{x+\frac{1}{2}{\mathrm{\Delta }}_x}=p_x+\frac{{\mathrm{\Delta }}_x}{2}f\left(q_x,p_{x+\frac{1}{2}{\mathrm{\Delta }}_x},u_x\right),(\mathrm{31a})$

$q_{x+{\mathrm{\Delta }}_x}=q_x+\frac{{\mathrm{\Delta }}_x}{2}\left[g\left(q_{x+{\mathrm{\Delta }}_x},p_{x+\frac{1}{2}\mathrm{\Delta }x},u_{x+\frac{1}{2}{\mathrm{\Delta }}_x}\right)+g\left(q_x,p_{x+\frac{1}{2}},u_{x+\frac{1}{2}{\mathrm{\Delta }}_x}\right)\right],(\mathrm{31b})$

$p_{x+{\mathrm{\Delta }}_x}=p_{x+\frac{1}{2}{\mathrm{\Delta }}_x}+\frac{{\mathrm{\Delta }}_x}{2}f\left(q_{x+{\mathrm{\Delta }}_x},p_{x+\frac{1}{2}{\mathrm{\Delta }}_x},u_{x+{\mathrm{\Delta }}_x}\right).(\mathrm{31c})$

The expression in Equations 29a, 29b can be retrieved by shifting Equation 31c back by ${\mathrm{\Delta }}_x$ and substituting the result into Equation 31a. Computationally, the form in Equations 29a, 29b is more efficient than Equations 31a, 31b, 31c (Hairer et al., 2003). If regular grid values are required, Equation 31c can be applied ad hoc. Unlike the non-general Störmer-Verlet scheme, the general scheme is not guaranteed to be explicit. If $f$ and $g$ are linear functions, then the resulting linearly implicit scheme can be resolved into an explicit update form.

When using the Störmer-Verlet method to discretize a PHS with a quadratic gradient, we no longer approximate the gradient. The output equation is untouched save for the addition of temporal indexes and averaging operators. Consider the quadratic PHS defined by a generally partitioned version of Equation 4

$\left[\begin{array}{c}\hfill {\partial }_t\mathbf{p}\hfill \\ \hfill {\partial }_t\mathit{q}\hfill \\ \hfill -\mathbf{y}\hfill \end{array}\right]=\left(\left[\begin{array}{ccc}\hfill {\mathbf{J}}_p\hfill & \hfill {\mathbf{K}}_{pq}\hfill & \hfill {\mathbf{G}}_p\hfill \\ \hfill -{\mathbf{K}}_{pq}^T\hfill & \hfill {\mathbf{J}}_q\hfill & \hfill {\mathbf{G}}_q\hfill \\ \hfill -{\mathbf{G}}_p^T\hfill & \hfill -{\mathbf{G}}_{\mathbf{q}}^T\hfill & \hfill -\mathbf{M}\hfill \end{array}\right]-\left[\begin{array}{ccc}\hfill {\mathbf{R}}_p\hfill & \hfill {\mathbf{R}}_{pq}\hfill & \hfill {\mathbf{P}}_p\hfill \\ \hfill {\mathbf{R}}_{qp}\hfill & \hfill {\mathbf{R}}_q\hfill & \hfill {\mathbf{P}}_q\hfill \\ \hfill {\mathbf{P}}_p^T\hfill & \hfill {\mathbf{P}}_{\mathbf{q}}^T\hfill & \hfill \mathbf{S}\hfill \end{array}\right]\right)\left[\begin{array}{c}\hfill {\mathbf{L}}_p\mathbf{p}\hfill \\ \hfill {\mathbf{L}}_q\mathit{q}\hfill \\ \hfill \mathbf{u}\hfill \end{array}\right],(32)$

where ${\mathbf{J}}_p$ and ${\mathbf{J}}_q$ are skew-symmetric matrices. All resistive matrices are $\ge 0$ and ${\mathbf{L}}_p$ and ${\mathbf{L}}_q$ are symmetric matrices arising from the quadratic Hamiltonian definition. A discretization with the Störmer-Verlet method produces the following system of equations. For brevity, we omit temporal indexes. $\mathbf{p}$ is located at times $n+\frac{1}{2}$ and $\mathit{q}$ and $\mathbf{u}$ are located at times $n$

${\delta }_{t-}\mathbf{p}=\left({\mathbf{J}}_p-{\mathbf{R}}_p\right)\left({\mathbf{L}}_p{\mu }_{t-}\mathbf{p}\right)+\left({\mathbf{K}}_{pq}-{\mathbf{R}}_{pq}\right)\left({\mathbf{L}}_q\mathit{q}\right)+\left({\mathbf{G}}_p-{\mathbf{P}}_p\right)\mathbf{u},(\mathrm{33a})$

${\delta }_{t+}\mathit{q}=\left({\mathbf{J}}_q-{\mathbf{R}}_q\right)\left({\mathbf{L}}_q{\mu }_{t+}\mathit{q}\right)+\left(-{\mathbf{K}}_{pq}^T-{\mathbf{R}}_{qp}\right)\left({\mathbf{L}}_p\mathbf{p}\right)+\left({\mathbf{G}}_q-{\mathbf{P}}_q\right){\mu }_{t+}\mathbf{u},(\mathrm{33b})$

$-\mathbf{y}=\left(-{\mathbf{G}}_p^T-{\mathbf{P}}_p^T\right)\left({\mu }_{t-}{\mathbf{L}}_p\mathbf{p}\right)+\left(-{\mathbf{G}}_q^T-{\mathbf{P}}_q^T\right)\left({\mathbf{L}}_q\mathit{q}\right)+\left(-\mathbf{M}-\mathbf{S}\right)\mathbf{u}.(\mathrm{33c})$

Any system discretized with the Störmer-Verlet method will be conditionally stable. We will not prove the general case and leave proofs for specific systems.

3.4 Energy and frequency analysis

Energy analysis is an informative tool that can aid the understanding of any FDTD scheme. In particular, it can help one determine if a scheme is conditionally or unconditionally stable, and—if a scheme is conditionally stable—what are the necessary criteria for stability. We evaluate the energy balance using the accumulated energy error, ${\mathfrak{h}}_{\mathit{err}}$, described in (Harrison-Harsley, 2018; Equation 3.75)

${\mathfrak{h}}_{\mathit{err}}^n=\frac{{\mathfrak{h}}^{n+1}-{\mathfrak{h}}^0+{\mathrm{\Delta }}_t{\sum }_{m=0}^n\left({\mathfrak{q}}^m-{\mathfrak{b}}^m\right)}{{\lfloor {\mathfrak{h}}^0\rfloor }_2},(34)$

where $\mathfrak{h}$, $\mathfrak{q}$ and $\mathfrak{b}$ are the discrete time correlates to the energetic quantities in Equation 5 and Equation 34 is a discrete-time integration of said equation. ${\lfloor \cdot \rfloor }_2$ denotes the rounding down to the nearest power of two used to reduce variations to machine precision. The stability of a scheme is established by deriving the conditions where ${\mathfrak{h}}^n,{\mathfrak{q}}^n\ge 0,\forall n$. In an unconditionally stable scheme this will always be true, but in conditionally stable schemes the limiting factors must be derived.

Frequency domain analysis is also a helpful tool and is equivalent to a $z$-Transform analysis (Bilbao, 2009). Frequency domain analysis is typically carried out using a frequency domain ansatz or single frequency solution $(u^n=z=e^{j\omega {\mathrm{\Delta }}_t})$ and is useful for evaluating the effects of numerical dispersion. Frequency domain analysis, in this manner, is limited to linear systems.

4 Example RLC circuit

An illustrative example of discrete-time PHS modeling is given by applying the framework towards modeling a series RLC circuit, shown in Figure 1. The RLC circuit appears throughout musical acoustics and is an equivalent model to a mechanical mass-spring-damper system which will be seen in the reed model (Section 5). Based on Kirchhoff’s voltage and current laws we have

$v_S=v_R+v_L+v_C,i_S=i_R=i_L=i_C,(35)$

and the system is governed by the following differential equation

$d_t{i}_s+2\alpha i_s+{\omega }_0^2\int i_{s}dt=\frac{1}{L}{v}_s,(36)$

Figure 1

Figure 1. Series RLC circuit. $q$ represents charge, $\varphi$ represents flux.

with $\alpha =\frac{R}{2L}$ representing damping in nepers and ${\omega }_0^2=\frac{1}{LC}$ the resonant frequency in rads/s. $d_t$ is the total derivative with respect to time. The Laplace domain admittance is

$Y\left(s\right)=\frac{I_s\left(s\right)}{V_s\left(s\right)}=\frac{\frac{1}{L}s}{s^2+2\alpha s+{\omega }_0^2}.(37)$

To form the PHS description, we first identify the power conjugate variables within the system. The variables corresponding to the storage elements are the voltage across and current through the inductor ($v_L$, $i_L$) and capacitor ($v_C$, $i_C$), the dissipative variables correspond to the voltage across and current through the resistor ($v_R$, $i_R$) and the external power variables are the voltage supplied by the ideal voltage source and the resulting current through the circuit. Based on Kirchhoff’s circuit laws and Ohm’s laws, the RLC system can be written as a PHS

$\left[\begin{array}{c}\hfill i_C=d_t\varphi \hfill \\ \hfill v_L=d_{t}q\hfill \end{array}\right]=\left((\right.\underset{\mathbf{J}}{\underbrace{\left[\begin{array}{cc}\hfill 0\hfill & \hfill -1\hfill \\ \hfill 1\hfill & \hfill 0\hfill \end{array}\right]}}-\underset{\mathbf{R}}{\underbrace{\left[\begin{array}{cc}\hfill R\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}}\left.\right))\left[\begin{array}{c}\hfill i_L=\frac{\partial }{\partial \varphi }\mathcal{H}\left(q,\varphi \right)\hfill \\ \hfill v_C=\frac{\partial }{\partial q}\mathcal{H}\left(q,\varphi \right)\hfill \end{array}\right]+\underset{\mathbf{G}}{\underbrace{\left[\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \end{array}\right]}}\left(v_s=u\right),(\mathrm{38a})$

$y=i_s=\underset{{\mathbf{G}}^T}{\underbrace{\left[\begin{array}{cc}\hfill 1\hfill & \hfill 0\hfill \end{array}\right]}}\left[\begin{array}{c}\hfill i_L=\frac{\partial }{\partial \varphi }\mathcal{H}\left(q,\varphi \right)\hfill \\ \hfill v_C=\frac{\partial }{\partial q}\mathcal{H}\left(q,\varphi \right)\hfill \end{array}\right].(\mathrm{38b})$

The internal state variables are $\mathbf{x}={\left[q\varphi \right]}^T$, the charge on the capacitor and the flux through the inductor. The total stored energy in the system is:

$\mathcal{H}\left(q,\varphi \right)=\frac{q^2}{2C}+\frac{{\varphi }^2}{2L},(39)$

and the Hamiltonian is quadratic. Notice that the elements describing a PHS from Equations 3a, 3b, $\mathbf{J}$, $\mathbf{R}$, and $\mathbf{G}$, arise naturally from physical laws.

4.1 Numerical schemes for the RLC circuit

We now apply the numerical methods discussed in Section 3 to the RLC PHS in Equation 4. Two criteria that are compared here are numerical stability and numerical dispersion. The accuracy of the two discretization methods is shown in Figure 2 in comparison to the analytical admittance given in Equation 37. Additionally, the accumulated energy error Equation 34 and a comparison of numerical dispersion in the two schemes is shown.

Figure 2

Figure 2. Top: admittance response of RLC PHS discretized with the discrete gradient (DG), discrete gradient with prewarping (DG-W), and Störmer-Verlet (SV) methods with a sampling rate of 48 kHz. The simulations are compared to the analytical response. Bottom: accumulated energy error for the DG scheme and SV scheme over the first 512 samples. Both display errors within machine precision, however, only the DG scheme conserves energy exactly. Right: Comparison of frequency warping of $f_0$ for the Störmer-Verlet and discrete gradient methods. All simulations were run for 4,096 samples.

4.1.1 Discrete gradient applied to RLC

The discrete gradient method is applied as described in Section 3. The discrete gradient RLC scheme is given by the following equations

${\delta }_{t+}{q}^n=\frac{1}{L}{\mu }_{t+}{\varphi }^n,{\delta }_{t+}{\varphi }^n=-\frac{1}{C}{\mu }_{t+}{q}^n-\frac{R}{L}{\mu }_{t+}{\varphi }^n-u^n,(\mathrm{40a})$

$y=-\frac{1}{L}{\mu }_{t+}{\varphi }^n.(\mathrm{40b})$

The scheme is unconditionally stable by definition of the discrete gradient method. We have the variation in discrete energy, $\mathfrak{h}$, is related to the energy at the boundary ports $\mathfrak{b}$ and dissipative ports $\mathfrak{q}$

${\delta }_{t+}\mathfrak{h}=\mathfrak{b}-\mathfrak{q},(41)$

with the discrete energies being an exact discretization of the continuous-time quantities,

$\mathfrak{h}=\frac{1}{2C}{\left(q^n\right)}^2+\frac{1}{2L}{\left({\varphi }^n\right)}^2,\mathfrak{b}=u^n{y}^n,\mathfrak{q}=R{\left({\mu }_{t+}\frac{{\varphi }^n}{L}\right)}^2.(42)$

For analyzing numerical dispersion, it is sufficient to analyze the unforced state-update equation, which leads to the characteristic equation

$z^2-2\frac{1-{\left(\frac{{\omega }_0{\mathrm{\Delta }}_t}{2}\right)}^2}{{\left(\frac{{\omega }_0{\mathrm{\Delta }}_t}{2}\right)}^2+\alpha {\mathrm{\Delta }}_t+1}z+\frac{{\left(\frac{{\omega }_0{\mathrm{\Delta }}_t}{2}\right)}^2-\alpha {\mathrm{\Delta }}_t+1}{{\left(\frac{{\omega }_0{\mathrm{\Delta }}_t}{2}\right)}^2+\alpha {\mathrm{\Delta }}_t+1}=0.(43)$

This can be compared to the analytical characteristic equation with $R_{\alpha }=e^{-\alpha {\mathrm{\Delta }}_t}$ and ${\omega }_n=\sqrt{{\alpha }^2-{\omega }_0^2}$

$z^2-2R_{\alpha }\cos \left({\omega }_n{\mathrm{\Delta }}_t\right)z+R_{\alpha }^2=0,(44)$

to evaluate the numerical dispersion of the scheme. For the RLC system, it is possible to adapt the parameters $L$ and $C$ such that the discretized system has the same characteristic equation as the analytical system (van Walstijn et al., 2016). This is a generalization of the “prewarping” procedure used when discretizing a system with the bilinear transform (Smith, 2010). The adjusted or prewarped parameters are

$\widehat{L}=\frac{R{\mathrm{\Delta }}_t}{4}\frac{1+2R_{\alpha }\cos \left({\omega }_n{\mathrm{\Delta }}_t\right)+R_{\alpha }^2}{1-R_{\alpha }^2},\widehat{C}=\frac{{\mathrm{\Delta }}_t}{R}\frac{1-R_{\alpha }^2}{1-2R_{\alpha }\cos \left({\omega }_n{\mathrm{\Delta }}_t\right)+R_{\alpha }^2}.(45)$

Figure 2 displays the response with and without prewarping. The discrete system response will deviate from the analytical response as ${\omega }_0$ increases in frequency without any prewarping.

4.1.2 Störmer-Verlet applied to RLC

The Störmer-Verlet method produces the following RLC scheme

${\delta }_{t-}{q}^{n+\frac{1}{2}}=\frac{1}{L}{\varphi }^n,{\delta }_{t+}{\varphi }^n=-\frac{1}{C}{q}^{n+\frac{1}{2}}-\frac{R}{L}{\mu }_{t+}{\varphi }^n-{\mu }_{t+}{u}^n,(\mathrm{46a})$

$y^n=-\frac{1}{L}{\varphi }^n.(\mathrm{46b})$

The scheme is conditionally stable, but not explicit. The exact discrete energy is not preserved and instead a modified energy is preserved:

${\delta }_{t+}{\mathfrak{h}}_{\text{mod}}=\mathfrak{b}-\mathfrak{q}.(47)$

Here,

${\mathfrak{h}}_{\text{mod}}=\frac{{\left({\varphi }^n\right)}^2}{2L}+\frac{{\left({\mu }_{t-}{q}^{n+\frac{1}{2}}\right)}^2}{2C}-\frac{{\mathrm{\Delta }}_t^2}{8L^{2}C}{\left({\varphi }^n\right)}^2,(48)$

and

$\mathfrak{q}=R{\left({\mu }_{t+}\frac{{\varphi }^n}{L}\right)}^2,\mathfrak{b}=\left({\mu }_{t+}{u}^n\right){\mu }_{t+}{y}^n.(49)$

The energetic quantities have been derived by taking the product of ${\mu }_{t+}\frac{{\varphi }^n}{L}$ with the second equation in Equation 46a and using the definitions in the first equation in Equation 46a, Equation 46b, and the identities in Equations 12, 16. The modified discrete energy is guaranteed to be positive if

$\mathfrak{h}\ge \frac{1}{2L}\left(1-\frac{{\mathrm{\Delta }}_t^2}{4LC}\right){\left({\varphi }^n\right)}^2.(50)$

For the discrete energy to remain positive ${\mathrm{\Delta }}_t\le 2\sqrt{LC}$ or, equivalently, $f_0\le \frac{f_s}{\pi }$ with $2\pi f_0={\omega }_0$. The characteristic equation for the unforced scheme is

$z^2+\frac{{\left({\omega }_0{\mathrm{\Delta }}_t\right)}^2-2}{1-\alpha {\mathrm{\Delta }}_t}z+\frac{1+\alpha {\mathrm{\Delta }}_t}{1-\alpha {\mathrm{\Delta }}_t}=0.(51)$

4.1.3 Comparison of different schemes

The RLC PHS was discretized with the discrete gradient method, with and without prewarping, and the Störmer-Verlet method. Results are shown in Figure 2. The system parameters are ${\omega }_0=2\pi 10$ kHz and $\alpha =0.15{\omega }_0$. The sampling rate is $f_s=48$ kHz. A high center frequency was chosen to test the accuracy of the schemes near the Nyquist limit. Both systems were used to process an impulse response 4,096 samples long. The stability of both schemes is demonstrated in the secondary plot comparing the accumulated energy error for the first 512 samples. There is no variation in the accumulated energy error following the initial transient response. Regarding accuracy, the discrete gradient method with prewarping (DG-W), performs the best. However, with no a priori of the system’s behavior, the Störmer-Verlet method (SV) performs slightly better than the discrete gradient method with no prewarping (DG). This result is further reinforced by comparing the frequency warping of $f_0$ for the discrete gradient and the Störmer-Verlet schemes using Equations 43, 51, respectively, over the range of stability for the Störmer-Verlet scheme. As shown on the right of Figure 2, generally, the Störmer-Verlet scheme has less frequency warping than the discrete gradient scheme below 12 kHz but the frequency warping increase dramatically as ${\omega }_0$ nears the limit of the stability criteria. Regarding memory and computational costs, we must store prior values of $u^n$ to compute ${\mu }_{t+}{u}^n$ and the average must also be computed. Both schemes require a linear system inverse and the general Störmer-Verlet scheme is not explicit.

5 The single reed model

The reed excitation mechanism acts as a blown closed pressure-controlled valve at the inlet of the woodwind bore (Fletcher, 1993). The reed model is nonlinear, an essential trait that produces steady oscillations from a continuous excitation (Chaigne and Kergomard, 2016). We will present a simple single-degree-of-freedom lumped model of the reed, a classical approach found in Backus (1963); Scavone (1997); van Walstijn and Avanzini (2007), among many others. The single reed has been modeled as a cantilever bar by Avanzini and van Walstijn (2004) and as a two-dimensional thin plate in Chatziioannou (2010). However, with a mind towards limiting computational complexity, only a lumped model is considered here.

A model of the single-reed consists of two nonlinear interactions: the reed contact with the lay called “beating”, a characteristic of single-reed instruments, and the nonlinear Bernoulli flow through the reed channel (Chaigne and Kergomard, 2016). We model the beating effect using the Hunt-Crossley impact model, previously utilized in Chatziioannou and van Walstijn (2012), and extend recent work on energy quadratization methods in musical acoustics (Bilbao et al., 2015b; Ducceschi et al., 2021; van Walstijn et al., 2024b) to develop a linearly implicit scheme. The system is modeled using the PHS framework. The PHS approach—not including contact dynamics, but including a turbulence model and 2D flow–was applied to modeling the lip-reed in brass instruments (Lopes and Hélie, 2016). In this section, we will first discuss the governing equations in the lumped reed model and then present a complete PHS description of the system. We will then derive discrete nonlinear and linear implicit schemes based on the discrete gradient method. Coupling the reed model to the bore model will be discussed later on in Section 8.

5.1 PHS model of the reed with contact dynamics

A lumped reed model aims to capture the behavior at the tip of the reed. We model the reed as a mass-spring-damper system which is a mechanical analogy of the RLC circuit presented in Section 4. The system is based on the position $y$ as seen in Figure 3. $y=0$ is the equilibrium position of the reed once the player has positioned their embouchure (the position and mechanical properties of the player’s lips), $y_l$ is the displacement from equilibrium to the mouthpiece lay, and $y_c$ is the displacement at which the Hunt-Crossley contact force begins to take effect. The lumped reed system is governed by the following equation:

$m{d}_t^{2}y+m\gamma d_{t}y+ky-f_c\left(h\right)=S_r{p}_{\mathrm{\Delta }},p_{\mathrm{\Delta }}=p_m-p_{in}.(52)$

$m$, $\gamma$, and $k$ are the effective mass, damping, and stiffness of the reed, which are affected by the player’s embouchure. $S_r$ is the effective reed area, also determined by the player’s embouchure, and $S_r{p}_{\mathrm{\Delta }}$ is the external force due to the pressure difference across the reed channel (van Walstijn and Avanzini, 2007). $p_m$ is the upstream pressure supplied by the player and $p_{in}$ is the pressure at the interface with the instrument bore. $f_c$ is the Hunt-Crossley contact force and is a function of compression: $h=y-y_c$. It is important to note that $y_c$ is not equivalent to $y_l$, as shown in Figure 3. $y_c$ has been empirically found to be a factor around $0.4-0.6$ times the displacement $y_l$ (Chatziioannou and van Walstijn, 2012; Chatziioannou et al., 2019). The contact force is defined as,

$f_c\left(h\right)=-{\partial }_h{\mathcal{V}}_c\left(h,t\right)-{\gamma }_c{\partial }_t{\mathcal{V}}_c\left(h,t\right)=-{\partial }_h{\mathcal{V}}_c\left(1+{\gamma }_c{\partial }_{t}h\right),(53)$

with ${\gamma }_c$ the contact damping and ${\mathcal{V}}_c$ the contact potential energy,

${\mathcal{V}}_c\left(h,t\right)=\frac{k_c}{\alpha +1}{\left[h\left(t\right)\right]}_{+}^{\alpha +1}.(54)$

$k_c$ represents the contact stiffness and $\alpha \ge 1$ is the contact power-law exponential. These parameters may be determined via the material parameters of the colliding objects (van Walstijn et al., 2024b) or through inverse modeling (Chatziioannou et al., 2019). The binary relationship between contact and no-contact is determined by the operator ${\left[\cdot \right]}_{+}$ where

${\left[h\right]}_{+}=\max \left(0,h\right).(55)$

Figure 3

Figure 3. Left: Diagram of the single reed and mouthpiece. Right: A lumped model of the reed including the Hunt-Crossley contact force.

The Hunt-Crossley model can be interpreted as a parallel combination of a nonlinear spring and damper that resists the incursion of one object into another.

A second nonlinearity within the reed model is the flow through the reed channel. We assume a steady, simple, and quasi-stationary Bernoulli flow model (Chaigne and Kergomard, 2016; Sec. 10.3.1.1)

$u_f=\text{sign}\left(p_{\mathrm{\Delta }}\right)S_j\left(y\right)\sqrt{\frac{2\left|p_{\mathrm{\Delta }}\right|}{\rho }},S_j\left(y\right)=w{\left[y_l-y\right]}_{+}.(56)$

$S_j(y)$ is the variable jet-area approximated by a rectangular surface with jet-width $w$. The steady Bernoulli flow model is found in Scavone (1997); van Walstijn and Avanzini (2007); Chatziioannou et al. (2019), among many others. It is likely pertinent to consider an unsteady flow model as was done by Lopes and Hélie (2016). However, we compensate for the simplified model by including a pumping flow,

$u_r=S_r{d}_{t}y,(57)$

that considers the contribution of reed motion to the flow. It was observed in Lopes and Hélie (2016) that the inclusion of a pumping flow is necessary for the model to be physically realistic and energy conserving. We now present a complete PHS description of the reed system based on Equations 52, 53, 56, 57, and denote momentum $\nu$, reserving $p$ for pressure variables

$\underset{{\mathbf{f}}_r}{\underbrace{d_t\left[\begin{array}{c}\hfill \nu \hfill \\ \hfill y\hfill \\ \hfill h\hfill \end{array}\right]}}=\left((\right.\underset{{\mathbf{J}}_r}{\underbrace{\left[\begin{array}{ccc}\hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}}-\underset{{\mathbf{R}}_r}{\underbrace{\left[\begin{array}{ccc}\hfill m\gamma +{\gamma }_c{\partial }_h{\mathcal{V}}_c\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}}\left.\right))\underset{{\mathbf{e}}_r}{\underbrace{\left[\begin{array}{c}\hfill {\partial }_{\nu }{\mathcal{T}}_r=d_{t}y\hfill \\ \hfill {\partial }_y{\mathcal{V}}_r=ky\hfill \\ \hfill {\partial }_h{\mathcal{V}}_c\hfill \end{array}\right]}}+\underset{{\mathbf{G}}_r}{\underbrace{\left[\begin{array}{cc}\hfill S_r\hfill & \hfill -S_r\hfill \\ \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}}\left[\begin{array}{c}\hfill p_m\hfill \\ \hfill p_{in}\hfill \end{array}\right],(\mathrm{58a})$

$\left[\begin{array}{c}\hfill u_m\hfill \\ \hfill u_{in}\hfill \end{array}\right]={\mathbf{G}}_r^T{\mathbf{e}}_r+\left[\begin{array}{c}\hfill u_f\hfill \\ \hfill -u_f\hfill \end{array}\right].(\mathrm{58b})$

$u_m$ and $u_{in}$ are the upstream and downstream volume velocities associated with $p_m$ and $p_{in}$, respectively. The total stored energy of the system in Equations 58a, 58b is

${\mathcal{H}}_r=\underset{{\mathcal{T}}_r}{\underbrace{\frac{{\nu }^2}{2m}}}+\underset{{\mathcal{V}}_r}{\underbrace{k\frac{y^2}{2}}}+{\mathcal{V}}_c.(59)$

It is worth noting that the flow out of the reed model $u_{in}=-(u_r+u_f)$ is negative in contrast to earlier literature due to sign conventions in defining the PHS. This is rectified when coupling the system with the woodwind bore. The pumping flow produced by ${\mathbf{G}}_r^T$ is a direct result of the PHS framework as the pumping flow complements the external forces due to $p_{\mathrm{\Delta }}$. It is not immediately clear that the system in Equation 58 maintains an energy balance due to the nonlinear contact damping and Bernoulli flow. Following the derivation in Bilbao et al. (2015b) and Chatziioannou et al. (2019), we can demonstrate passivity by taking the inner product of the flows and efforts in Equation 58a

$⟨{\mathbf{e}}_r,{\mathbf{f}}_r⟩={\partial }_t{\mathcal{H}}_r=-\underset{{\mathcal{Q}}_r}{\underbrace{m\gamma {\left({\partial }_{\nu }\mathcal{T}\right)}^2}}-\underset{{\mathcal{Q}}_c}{\underbrace{\left({\gamma }_c{\partial }_h{\mathcal{V}}_c\right){\left({\partial }_{\nu }\mathcal{T}\right)}^2}}+⟨{\mathbf{e}}_r,{\mathbf{G}}_r{\left[p_m{p}_{in}\right]}^T⟩,(60)$

and identifying ${\mathcal{Q}}_r$ as the damping due to the lumped reed model and ${\mathcal{Q}}_c$ as the damping due to the contact force. ${\mathcal{Q}}_c$ is non-negative as ${\partial }_h{\mathcal{V}}_c$ is non-negative by definition. To simplify the last product in Equation 60, the inner product of Equation 58b is taken with ${\left[p_m{p}_{in}\right]}^T$ resulting in the following energetic quantities

$\begin{array}{ll}\hfill ⟨{\mathbf{e}}_r,{\mathbf{G}}_r{\left[p_m{p}_{in}\right]}^T⟩& =p_{in}{u}_{in}+p_m{u}_m-p_m{u}_f+p_{in}{u}_f,\hfill \\ \hfill & =\underset{{\mathcal{B}}_r}{\underbrace{p_{in}{u}_{in}+p_m{u}_m}}-\underset{{\mathcal{Q}}_f}{\underbrace{u_f{p}_{\mathrm{\Delta }}}}.\hfill \end{array}(61)$

${\mathcal{B}}_r$ is the power flowing through the boundaries and ${\mathcal{Q}}_f$ is the dissipation due to the Bernoulli flow through the reed channel equal to

${\mathcal{Q}}_f=S_j\text{sign}\left(p_{\mathrm{\Delta }}\right)\sqrt{\frac{2}{\rho }}{\left|p_{\mathrm{\Delta }}\right|}^{1/2}\cdot p_{\mathrm{\Delta }}=S_j\sqrt{\frac{2}{\rho }}{\left|p_{\mathrm{\Delta }}\right|}^{3/2}\ge 0.(62)$

The dissipation due to the Bernoulli flow is guaranteed to be non-negative by $p_{\mathrm{\Delta }}=\text{sign}\left(p_{\mathrm{\Delta }}\right)\left|p_{\mathrm{\Delta }}\right|$ and $\text{sign}{\left(x\right)}^2=1$. Restating the energy balance of the system,

${\partial }_t{\mathcal{H}}_r={\mathcal{B}}_r-\underset{\ge 0}{\underbrace{\left({\mathcal{Q}}_r+{\mathcal{Q}}_c+{\mathcal{Q}}_f\right)}},(63)$

the system in Equations 58a, 58b maintains an energy balance as all the dissipative terms and ${\mathcal{H}}_r$ are all guaranteed to be non-negative.

5.2 A nonlinear implicit scheme

We discretize the system in Equations 58a, 58b using the discrete gradient method and obtain Equation 64 by approximation of the time derivative and gradients. This results in a pair of implicit nonlinear equations in both the state-update equation (due to the Hunt-Crossley force) and the output equation (due to the Bernoulli flow). We will focus on the state-update equation and the discretization of the output equation will be left to Section 8. We are particularly concerned with the discrete-time transcendental equation:

${\delta }_{t+}{\nu }^n=-k{\mu }_{t+}{y}^n-{\delta }_{h+}{\mathcal{V}}_c^n\left(1+{\gamma }_c{\mu }_{t+}\left(\frac{{\nu }^n}{m}\right)\right)-\gamma m{\mu }_{t+}\left(\frac{{\nu }^n}{m}\right)+S_r{p}_{\mathrm{\Delta }},(64)$

which corresponds to the discretized form of the first line in Equation 58a. It is helpful here to introduce an intermediary variable:

$\xi =h^{n+1}-h^n=y^{n+1}-y^n={\mathrm{\Delta }}_t{\delta }_{t+}y.(65)$

We apply the property in Equation 13 to replace ${\delta }_{t+}{\nu }^n$ and ${\mu }_{t+}{y}^n$ with factors of ${\mu }_{t+}{\nu }^n$ and ${\delta }_{t+}{y}^n$, respectively. Then, we utilize the discretized form of the second line in Equation 58a,

${\delta }_{t+}{y}^n={\mu }_{t+}\frac{{\nu }^n}{m},(66)$

and Equation 65 to rewrite the transcendental equation as a function whose root is equal to $\xi$

$F\left(\xi \right)=\xi -{\xi }^{\ast }+\frac{1}{a_0}{\delta }_{h+}{\mathcal{V}}_c^n\left(1+\frac{{\gamma }_c}{{\mathrm{\Delta }}_t}\xi \right)=0,(67)$

with

${\xi }^{\ast }=\frac{1}{a_0}\left(S_r{p}_{\mathrm{\Delta }}+a_1^n\right),(\mathrm{68a})$

$a_0=\frac{2m}{{\mathrm{\Delta }}_t^2}+\frac{k}{2}+\frac{m\gamma }{{\mathrm{\Delta }}_t},(\mathrm{68b})$

$a_1^n=\frac{2}{{\mathrm{\Delta }}_t}{\nu }^n-k{y}^n.(\mathrm{68c})$

${\xi }^{\ast }$ represents the solution to $\xi$ in the absence of the contact force. We approximate ${\delta }_{h+}{\mathcal{V}}_c$ using the definition of the discrete gradient in Equation 19,

$F\left(\xi \right)=\xi -{\xi }^{\ast }+\frac{1}{a_0}\frac{{\mathcal{V}}_c\left(\xi +h^n\right)-{\mathcal{V}}_c\left(h^n\right)}{\xi }\left(1+\frac{{\gamma }_c}{{\mathrm{\Delta }}_t}\xi \right).(69)$

This equation can be solved using a Newton-Raphson iteration with a guaranteed solution due to the convexity of $F(\xi )$ (Bilbao et al., 2015b). Division by $\xi$ can prevent convergence and a more numerically robust definition of $F(\xi )$ and ${\partial }_{\xi }F(\xi )$ can be produced following the procedure in van Walstijn et al. (2024b). After $\xi$ has been derived, the next states ${\nu }^{n+1}$, $y^{n+1}$, and $h^{n+1}$ can be updated based on $\xi$ via

${\nu }^{n+1}=\frac{2m}{{\mathrm{\Delta }}_t}\xi -{\nu }^n,(\mathrm{70a})$

$y^{n+1}=\xi +y^n,(\mathrm{70b})$

$h^{n+1}=\xi +h^n.(\mathrm{70c})$

The discrete Hamiltonian and dissipated energy terms associated with this scheme are

${\mathfrak{h}}_r=\frac{{\left({\nu }^n\right)}^2}{2m}+k\frac{{\left(y^n\right)}^2}{2}+{\mathfrak{h}}_c,{\mathfrak{h}}_c={\mathcal{V}}_c\left(h^n\right),(\mathrm{71a})$

${\mathfrak{h}}_r=\left(m\gamma +{\gamma }_c{\partial }_h{\mathcal{V}}_c\left(h^n\right)\right){\left({\mu }_{t+}\frac{{\nu }^n}{m}\right)}^2.(\mathrm{71b})$

5.3 A linearly implicit scheme

Explicit and linearly implicit schemes representing Hertzian contact dynamics have been derived using energy quadratization strategies in Ducceschi et al. (2021) and van Walstijn et al. (2024b). Explicit Hunt-Crossley contact dynamics were recently presented in van Walstijn et al. (2024a). Energy quadratization has also been used in the discretization of nonlinear port-Hamiltonian systems (Lopes et al., 2015). Invariant energy quadratization (IEQ) (Zhao et al., 2016) and scalar auxiliary variable (SAV) methods (Shen et al., 2018) apply a change of basis to the energy function by representing the energy with a quadratic variable. For single-variable problems, as treated here, the methods are equivalent. We will only quadratize the contact potential, leaving the other energies untouched. The contact potential is now defined based on an auxiliary state $\sigma$ with energy ${\mathcal{V}}_{\sigma }$,

${\mathcal{V}}_{\sigma }=\frac{1}{2}{\sigma }^2,(72)$

and $h$ is related to the new state by the function $f(\cdot )$ with $\sigma =f(h)$. This substitution is viable here as the given energy function is non-negative (Shen et al., 2018). ${\partial }_h{\mathcal{V}}_c$ and ${\partial }_{t}h$ are defined relative to the auxillary state by,

${\partial }_h{\mathcal{V}}_c=g\sigma ,{\partial }_t\sigma =g{\partial }_{t}h,(73)$

where $g$ is the gradient variable equal to ${\partial }_{h}f(h)$,

$g=\sqrt{\frac{1}{2}{k}_c\left(\alpha +1\right){\left[h\right]}_{+}^{\alpha -1}}.(74)$

Applying Equation 73 to the PHS in Equations 58a, 58b we derive a new PHS with elements

${\mathbf{J}}_{\sigma }=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill -1\hfill & \hfill -g\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill g\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right],{\mathbf{R}}_{\sigma }=\left[\begin{array}{ccc}\hfill m\gamma +{\gamma }_{c}g\sigma \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right],(\mathrm{75a})$

and efforts and flows,

${\mathbf{e}}_{\sigma }=\left[\begin{array}{c}\hfill {\partial }_{\nu }{\mathcal{T}}_r\hfill \\ \hfill {\partial }_y{\mathcal{V}}_r\hfill \\ \hfill {\partial }_{\sigma }{\mathcal{V}}_{\sigma }\hfill \end{array}\right],{\mathbf{f}}_{\sigma }=d_t\left[\begin{array}{c}\hfill \nu \hfill \\ \hfill y\hfill \\ \hfill \sigma \hfill \end{array}\right].(\mathrm{75b})$

The total stored energy in the system is now

${\mathcal{H}}_{\sigma }={\mathcal{T}}_r+{\mathcal{V}}_r+{\mathcal{V}}_{\sigma }.(\mathrm{75c})$

Following Ducceschi et al. (2021), we discretize the system by placing ${\sigma }^{n+\frac{1}{2}}$ on the interleaved $n+\frac{1}{2}$ time step such that ${\sigma }^{n+\frac{1}{2}}$ is an independent time series. As a result of the interleaved step, we have an increase in memory requirements as we must store the states ${\sigma }^{n+\frac{1}{2}}$ and ${\sigma }^{n-\frac{1}{2}}$. We use the following approximations in our discretization

${\partial }_{\sigma }{\mathcal{V}}_{\sigma }\approx {\mu }_{t\cdot }{\sigma }^{n+\frac{1}{2}},{\partial }_h{\mathcal{V}}_c=g\sigma \approx g^n{\sigma }^{n-\frac{1}{2}},(\mathrm{76a})$

${\partial }_t\sigma \approx {\delta }_{t\cdot }{\sigma }^{n+\frac{1}{2}}=g^n{\delta }_{t+}{y}^n,(\mathrm{76b})$

where $g^n$ is a discretization of Equation 74. Using Equation 76a, the discrete form of Equation 67 is rewritten as

$F\left(\xi \right)=\xi -{\xi }^{\ast }+\frac{1}{a_0}{g}^n{\mu }_{t\cdot }{\sigma }^{n+\frac{1}{2}}+{\sigma }^{n-\frac{1}{2}}\left(\frac{{\gamma }_c}{{\mathrm{\Delta }}_t}\xi \right)=0,(77)$

and $\xi$ can be solved using the property in Equation 14 and Equation 76b

$\xi =\frac{{\xi }^{\ast }-\frac{g^n}{a_0}{\sigma }^{n-\frac{1}{2}}}{1+\frac{g^n}{a_0}\left(g^n+\frac{{\gamma }_c}{{\mathrm{\Delta }}_t}{\sigma }^{n-\frac{1}{2}}\right)}.(78)$

As discussed by Ducceschi et al. (2021) and van Walstijn et al. (2024b) the value of $g^n$ must be chosen so that the sign of the force is consistent and artifacts related to the contact energy are non-zero in the absence of contact. The former issue is handled by setting $g^n$ to have the same sign as ${\sigma }^{n-\frac{1}{2}}$. This also ensures that the contact damping coefficient, ${\gamma }_c{g}^n{\sigma }^{n-\frac{1}{2}}$, is always positive. The latter issue, pertaining to non-zero contact energy in the absence of contact, is derived via the discretized form of the second equation in Equation 73 for ${\sigma }^{n+\frac{3}{2}}=0$. As a result, $g^n$ is defined as

$g^n=\left\{\begin{array}{cc}-\frac{1}{2}\frac{{\sigma }^{n-\frac{1}{2}}}{{\xi }^{\ast }},\hfill & h^n<0\hfill \\ \text{sign}\left({\sigma }^{n-\frac{1}{2}}\right)\sqrt{\frac{1}{2}{k}_c\left(\alpha +1\right){\left[h^n\right]}_{+}^{\alpha -1}},\hfill & h^n\ge 0\hfill \end{array}\right..(79)$

Following van Walstijn et al. (2024b), a third branch would be given by determining under what conditions, in Equation 78, $\xi ={\xi }^{\ast }$. However, in the authors’ experiments it was found that including this condition produces an undesirable transient when $\xi$ is preemptively set to ${\xi }^{\ast }$ and so a third branch is omitted from the model presented here. Finally, $\sigma$ can be updated using Equation 76b

${\sigma }^{n+\frac{3}{2}}={\sigma }^{n-\frac{1}{2}}+2g^n\xi .(80)$

The discrete variation in energy associated with the linearly implicit scheme is derived using the properties of Equations 12, 15

${\delta }_{t+}{\mathfrak{h}}_{\sigma }^n={\delta }_{t+}\left[{\mu }_{t-}\left(\frac{{\left({\sigma }^{n+\frac{1}{2}}\right)}^2}{2}\right)\right],(\mathrm{81a})$

with $\sigma$ replacing ${\mathfrak{h}}_c^n$ in Equation 71a. The discrete dissipated energy is now

${\mathfrak{q}}_{\sigma }^n=\left(m\gamma +{\gamma }_c{g}^n{\sigma }^{n-\frac{1}{2}}\right){\left({\mu }_{t+}\frac{{\nu }^n}{m}\right)}^2.(\mathrm{81b})$

In Figure 4, we compare the hysteretic compressive behavior associated with the Hunt-Crossley contact force for the nonlinear and linearly implicit schemes with different initial velocities. As can be seen, the linearly implicit scheme follows the behavior of the nonlinear implicit scheme but does not match it exactly due to how the energy quadratized PHS is discretized. The accuracy decreases with increasing velocity.

Figure 4

Figure 4. Hunt-Crossley force versus deformation simulated for different initial velocities 0.5, 1, and 2 m/s using the nonlinear implicit scheme (blue) and the linearly implicit scheme (red, dashed). Parameters used for the plot are $m=100$ mg, $k=10$ N/m, $\gamma =100{\text{s}}^{-1}$, $k_c=1\times 10^7{\text{N/m}}^{\alpha }$, ${\gamma }_c=5$ s/m, and $\alpha =2.5$.

6 Woodwind bore model

The woodwind bore can be modeled in one-dimension as a pipe with a varying cross-sectional area. Viscothermal losses are incorporated according to the model of Zwikker and Kosten (1949). We represent these losses in our model using the network synthesis method proposed in Bilbao et al. (2015a). A PHS formulation of the Webster-Lokshin model which approximates viscothermal losses with various simplifications was proposed in Gorrec and Matignon (2013). We will summarize the frequency domain description of the lossy horn equation, then present a PHS representation of a network synthesis approximation of the viscothermal losses and the overall system. To reduce frequency warping, we will discretize the system with the Störmer-Verlet method. In the absence of losses, the resulting scheme is exactly the scheme presented in Bilbao and Harrison (2016) though derived from a different perspective.

6.1 Frequency-domain description of the horn equation

We consider a tube of varying circular cross-sectional area of length $L$ defined on the domain $\mathcal{L}=0\le z\le L$ where $z$ is the symmetric axis down the length of the bore. A frequency domain model including the effects of viscothermal losses is (Caussé et al., 1984; Keefe, 1984)

${\partial }_{z}P+j\omega \frac{\rho }{S\left(z\right)}U+Z_{v}U=0,(\mathrm{82a})$

${\partial }_{z}U+j\omega \frac{S\left(z\right)}{\rho c^2}P+Y_{t}P=0.(\mathrm{82b})$

Here, $P$ and $U$ are the Fourier transforms of pressure and volume velocity, respectively. These equations are the first-order form of the horn equation (Chaigne and Kergomard, 2016). Instead of modeling $P$ and $U$ directly, we will apply a change of variables,

$\mathrm{\Psi }=\frac{S\left(z\right)}{\rho c^2}P,\mathrm{\Upsilon }=\frac{\rho }{S\left(z\right)}U,(83)$

producing the equivalent system

${\partial }_z\left(\frac{\rho c^2}{S\left(z\right)}\mathrm{\Psi }\right)+j\omega \mathrm{\Upsilon }+Z_v\frac{S\left(z\right)}{\rho }\mathrm{\Upsilon }=0,(\mathrm{84a})$

${\partial }_z\left(\frac{S\left(z\right)}{\rho }\mathrm{\Upsilon }\right)+j\omega \mathrm{\Psi }+Y_t\frac{\rho c^2}{S\left(z\right)}\mathrm{\Psi }=0.(\mathrm{84b})$

$\mathrm{\Psi }$ and $\mathrm{\Upsilon }$ are the underlying state variables related to pressure and volume velocity, respectively.

$Z_v$ and $Y_t$ are frequency dependent functions encapsulating viscous and thermal losses and are defined as (Caussé et al., 1984; Bilbao and Harrison, 2016)

$Z_v=j\omega \frac{\rho }{S\left(z\right)}\frac{F_v}{1-F_v},Y_t=j\omega \frac{S\left(z\right)}{\rho c^2}\left(\gamma -1\right)F_t,(85)$

with

$F_v=2\frac{J_1\left({\eta }_v\right)}{{\eta }_v{J}_0\left({\eta }_v\right)},F_t=2\frac{J_1\left({\eta }_t\right)}{{\eta }_t{J}_0\left({\eta }_t\right)}(\mathrm{86a})$

${\eta }_v=r\left(z\right)\sqrt{\frac{j\omega }{c{l}_v}},l_v=\frac{\mu }{\rho c},(\mathrm{86b})$

${\eta }_t=r\left(z\right)\sqrt{\frac{j\omega }{c{l}_t}},l_t=\frac{l_v}{P_r}.(\mathrm{86c})$

$\gamma$ is the ratio of specific heats. $J_0\left(\cdot \right)$ and $J_1\left(\cdot \right)$ are zeroth and first-order Bessel functions of the first kind. $r(z)$ is the radius at each position $z$, $c$ is the speed of sound, and $l_v$ and $l_t$ are the viscous and thermal boundary layer thicknesses. $\mu$ is first coefficient of viscosity of air and $\rho$ is the density of air. $P_r=\frac{C_p\mu }{\kappa }$ is the Prandtl number with $C_p$ the coefficient of specific heat at constant pressure and $\kappa$ is the coefficient of thermal conductivity. Keefe (1984) provides equations for calculating the value of these constants within $\pm 10*$C of $26.85*$C. A time-domain form of Equation 84a is

${\partial }_z\left(\frac{\rho c^2}{S\left(z\right)}\psi \right)+{\partial }_t\upsilon +\underset{p_v}{\underbrace{z_v\ast \left(\frac{S\left(z\right)}{\rho }\upsilon \right)}}=0,(\mathrm{87a})$

${\partial }_z\left(\frac{S\left(z\right)}{\rho }\upsilon \right)+{\partial }_t\psi +\underset{u_t}{\underbrace{y_t\ast \left(\frac{\rho c^2}{S\left(z\right)}\psi \right)}}=0,(\mathrm{87b})$

where $\ast$ operator denotes convolution. $p_v$ and $u_t$ are variables representing the effects of viscous and thermal losses, respectively. $p_v$ and $u_t$ can be approximated using network synthesis (Bilbao et al., 2015a). We will interpret $p_v$ and $u_t$ as the output of PHS loss sub-models.

6.2 Viscothermal loss approximation

$Z_v$ and $Y_t$ are approximated through an $M$ section Foster network synthesis with the form shown in Figure 5. The frequency-domain transfer function for the losses at each spatial location are

$Z_v=R_0+{\displaystyle \sum _{m=1}^M}\frac{R_{m}j\omega }{j\omega +R_m/L_m},(\mathrm{88a})$

$Y_t={\displaystyle \sum _{m=1}^M}\frac{1/G_{m}j\omega }{j\omega +1/\left(G_m{C}_m\right)}.(\mathrm{88b})$

Figure 5

Figure 5. $M$-branch RL Foster I structure approximating $Z_v$ (top) and $M$-branch RC Foster II structure approximating $Y_t$ (bottom).

$R_0$ is the limit as $Z_v$ approaches zero frequency. $Y_t(0)=0$ and no constant term is needed (Chatziioannou et al., 2019). Both of these systems can be represented as $M^{\text{th}}$-order PHSs. The losses due to viscous effects, $p_v$, are approximated by the following PHS

${\nabla }_{\mathit{\varphi }}{\mathcal{H}}_v=-{\mathbf{R}}_v^{-1}{\partial }_t\mathit{\varphi }+\mathbf{G}\left(\frac{S\left(z\right)}{\rho }\upsilon \right),(\mathrm{89a})$

$p_v={\mathbf{G}}^T{\partial }_t\mathit{\varphi }+R_0\left(\frac{S\left(z\right)}{\rho }\upsilon \right),(\mathrm{89b})$

with fluxes $\mathit{\varphi }={\left[{\varphi }_1,\dots ,{\varphi }_M\right]}^T$. ${\mathbf{R}}_v$ is an $M\times M$ diagonal resistance matrix and $\mathbf{G}$ is an $M\times 1$ vector of ones. These elements are defined as:

${\mathbf{R}}_v=\text{diag}\left(\left[R_1,\dots ,R_M\right]\right),\mathbf{G}={\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \dots \hfill & \hfill 1\hfill \end{array}\right]}^T,{\mathcal{H}}_v\left({\varphi }_m\right)=\frac{{\varphi }_m^2}{2L_m},\forall {\varphi }_m\in \mathit{\varphi }.(\mathrm{89c})$

The losses due to thermal effects, $u_t$, are approximated by the PHS

${\nabla }_{\mathit{q}}{\mathcal{H}}_t=-{\mathbf{R}}_t{\partial }_t\mathit{q}+\mathbf{G}\left(\frac{\rho c^2}{S\left(z\right)}\psi \right),(\mathrm{90a})$

$u_t={\mathbf{G}}^T{\partial }_t\mathit{q},(\mathrm{90b})$

with charges $\mathit{q}={\left[q_1,\dots ,q_M\right]}^T$. ${\mathbf{R}}_t$ is an $M\times M$ diagonal resistance matrix and $\mathbf{G}$ is an $M\times 1$ vector of ones. These elements are defined as:

${\mathbf{R}}_t=\text{diag}\left(\left[1/G_1,\dots ,1/G_M\right]\right),\mathbf{G}={\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \dots \hfill & \hfill 1\hfill \end{array}\right]}^T,{\mathcal{H}}_t\left(q_m\right)=\frac{q_m^2}{2C_m},\forall {\varphi }_m\in \mathit{q}.(\mathrm{90c})$

Both systems are defined as conjugate PHS systems as the output of both systems are the sum of the flows of their respective storage components. As described in Section 3.2, discretization with the discrete gradient method results in a linearly implicit scheme.

6.3 Power balance of the horn model including loss approximation

Following Bilbao and Harrison (2016), a power balance of the horn model, including the loss approximations, is derived by taking the product $\frac{S(z)}{\rho }\upsilon$ with Equation 87a and $\frac{\rho c^2}{S(z)}\psi$ with Equation 87b. Integrating over the domain $\mathcal{L}$ and summing the equations,

${\partial }_t\left({\int }_{\mathcal{L}}{\mathcal{H}}_{h}dz\right)-\mathcal{B}+{\int }_{\mathcal{L}}{p}_v\frac{S\left(z\right)}{\rho }\upsilon dz+{\int }_{\mathcal{L}}{u}_t\frac{\rho c^2}{S\left(z\right)}\psi dz=0.(91)$

${\mathcal{H}}_h$ represents the energy storage density in absence of losses and $\mathcal{B}$ is the power supplied at the bore boundaries,

${\mathcal{H}}_h=\frac{1}{2}\frac{\rho c^2}{S\left(z\right)}{\psi }^2+\frac{1}{2}\frac{S\left(z\right)}{\rho }{\upsilon }^2,\mathcal{B}=p\left(0\right)u\left(0\right)-p\left(L\right)u\left(L\right).(92)$

We can identify the quantities $L_{\psi }(z)=\frac{\rho c^2}{S(z)}$ and $L_{\upsilon }(z)=\frac{S(z)}{\rho }$ related to the quadratic Hamilton in Equation 92. The terms related to viscothermal losses in Equation 91 are equivalent to the product of the external ports of the PHS sub-models in Equations 89a, 89b and Equations 90a, 90b, owing to their definition as PHSs. As such, we can directly use the power balance property from Equation 5 to simplify our derivation of the power balance. The power balance of the entire horn system including losses is

${\partial }_t\left[{\int }_{\mathcal{L}}\left({\mathcal{H}}_h+{\mathcal{H}}_v+{\mathcal{H}}_t\right)dz\right]=\mathcal{B}-{\int }_{\mathcal{L}}\left({\mathcal{Q}}_v+{\mathcal{Q}}_t\right)dz,(93)$

with

${\mathcal{Q}}_v=⟨{\partial }_t\mathit{\varphi },{\mathbf{R}}_v^{-1}{\partial }_t\mathit{\varphi }⟩+R_0\Vert L_{\upsilon }\left(z\right)\upsilon {\Vert }^2,{\mathcal{Q}}_t=⟨{\partial }_t\mathit{q},{\mathbf{R}}_t{\partial }_t\mathit{q}⟩.(94)$

The total system can be written as an infinite dimensional PHS based on the variational derivatives:

$\frac{\delta {\mathcal{H}}_h}{\delta \psi }={\partial }_{\psi }{\mathcal{H}}_h=L_{\psi }\left(z\right)\psi =p,\frac{\delta {\mathcal{H}}_h}{\delta \upsilon }={\partial }_{\upsilon }{\mathcal{H}}_h=L_{\upsilon }\left(z\right)\upsilon =u.(95)$

The variational derivatives are equivalent to the derivatives of the energy density when the energy density does not depend on the derivatives of the state (Duindam et al., 2009). Here, we use Leibniz’s notation in conjunction with the operator $\delta$ to refer specifically to the variational derivative. To prevent confusion, we will use ${\partial }_{\psi }{\mathcal{H}}_h$ and ${\partial }_{\upsilon }{\mathcal{H}}_h$ going forward. Due to the change in variables in Equation 83, the variational derivatives are equal to the pressure and volume velocity variables at each spatial location. The main benefit of this change in variables is that the boundary efforts and flows in Equation 94 correspond directly with pressure and volume velocity which are the acoustic domain power-conjugated variables.

The distributed PHS formulation of the horn equation with viscothermal losses is then:

$\left[\begin{array}{c}\hfill {\partial }_t\upsilon \hfill \\ \hfill {\partial }_t\psi \hfill \\ \hfill {\nabla }_{\mathit{\varphi }}{\mathcal{H}}_v\hfill \\ \hfill {\nabla }_{\mathit{q}}{\mathcal{H}}_t\hfill \end{array}\right]=\left(\underset{\mathcal{J}}{\underbrace{\left[\begin{array}{cccc}\hfill 0\hfill & \hfill -{\partial }_z\hfill & \hfill -{\mathbf{G}}^T\hfill & \hfill 0\hfill \\ \hfill -{\partial }_z\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{\mathbf{G}}^T\hfill \\ \hfill \mathbf{G}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathbf{G}\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}}\right.-\left.\underset{\mathcal{R}}{\underbrace{\left[\begin{array}{cccc}\hfill R_0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathbf{R}}_v\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathbf{R}}_t\hfill \end{array}\right]}}\right)\left[\begin{array}{c}\hfill {\partial }_{\upsilon }{\mathcal{H}}_h\hfill \\ \hfill {\partial }_{\psi }{\mathcal{H}}_h\hfill \\ \hfill {\partial }_t\mathit{\varphi }\hfill \\ \hfill {\partial }_t\mathit{q}\hfill \end{array}\right].(96)$

The upper left quadrant in $\mathcal{J}$ — representing the interconnection of efforts and flows in the lossless horn equation—is skew-symmetric by the inclusion of the boundary elements $p(0)$, $u(0)$, $p(L)$, and $u(L)$ (Duindam et al., 2009).

6.4 Discretization of the horn equation by the Störmer-Verlet method

In light of the frequency warping associated with the discrete gradient method we now discretize the system in Equation 96 in time and space with the Störmer-Verlet method. This produces a scheme that, excluding viscothermal losses, is equivalent to the scheme proposed in Bilbao and Harrison (2016). Spatial discretization with the Störmer-Verlet scheme maintains the skew-symmetry property of the upper left block matrix in Equation 96. However, this property combined with the reduced time discretization error does completely characterize the accuracy of the scheme. Analysis of the scheme accuracy in relation to the symplectic property of the discretization method merits further investigation.

The fully discretized horn equation scheme is

${\delta }_{t-}{\upsilon }_{l+\frac{1}{2}}^{n+\frac{1}{2}}+{\delta }_{z+}\left(\frac{\rho c^2}{S_l}{\psi }_l^n\right)+p_{v,l+\frac{1}{2}}^n=0,(\mathrm{97a})$

${\delta }_{t+}{\psi }_l^n+{\delta }_{z-}\left(\frac{S_{l+\frac{1}{2}}}{\rho }{\upsilon }_{l+\frac{1}{2}}^{n+\frac{1}{2}}\right)+u_{t,l}^{n+\frac{1}{2}}=0.(\mathrm{97b})$

Importantly, we define $S_l={\mu }_{z-}{S}_{l+\frac{1}{2}}$ and the bore is sampled on the interleaved grid based on an energetic stability analysis of the scheme (Harrison-Harsley, 2018). Equation 97b is centered on the interleaved grid in $n+\frac{1}{2}$ time and the regular grid $l$ in space. Equation 97a is defined in the opposite manner. The output of the losses are then defined accordingly. The final step of the Störmer-Verlet method in Equation 31c allows us to compute the values of $\upsilon$ at boundaries $l=0$ and $l=N$

$\frac{S_{\frac{1}{2}}}{\rho }{\upsilon }_{\frac{1}{2}}^{n+\frac{1}{2}}=u_0^{n+\frac{1}{2}}-\frac{{\mathrm{\Delta }}_z}{2}\left[{\delta }_{t+}{\psi }_0^n+u_{t,0}^{n+\frac{1}{2}}\right],(\mathrm{98a})$

$u_N^{n+\frac{1}{2}}=\frac{S_{N-\frac{1}{2}}}{\rho }{\upsilon }_{N-\frac{1}{2}}^{n+\frac{1}{2}}-\frac{{\mathrm{\Delta }}_z}{2}\left[{\delta }_{t+}{\psi }_N^n+u_{t,N}^{n+\frac{1}{2}}\right].(\mathrm{98b})$

This equivalence can also be derived by substituting the property ${\delta }_{z-}=\frac{2}{{\mathrm{\Delta }}_z}\left({\mu }_{z-}-e_{z-}\right)$ into Equation 97b at $l=0$ and $l=N$. The viscothermal losses are discretized using the discrete gradient method at each spatial step

${\mu }_{t-}\left({\mathbf{L}}_v{\mathit{\varphi }}^n\right)={\mathbf{R}}_v^{-1}{\delta }_{t-}{\mathit{\varphi }}^n+\mathbf{G}{\mu }_{t-}\left(\frac{S}{\rho }{\upsilon }^{n+\frac{1}{2}}\right),(\mathrm{99a})$

$p_v^n={\mathbf{G}}^T{\delta }_{t-}{\mathit{\varphi }}^n+R_0{\mu }_{t-}\left(\frac{S}{\rho }{\upsilon }^{n+\frac{1}{2}}\right),(\mathrm{99b})$

${\mu }_{t+}\left({\mathbf{L}}_t{\mathit{q}}^{n+\frac{1}{2}}\right)={\mathbf{R}}_t{\delta }_{t+}{\mathit{q}}^{n+\frac{1}{2}}+\mathbf{G}{\mu }_{t+}\left(\frac{\rho c^2}{S}{\psi }^n\right),(\mathrm{100a})$

$u_t^{n+\frac{1}{2}}={\mathbf{G}}^T{\delta }_{t+}{\mathit{q}}^{n+\frac{1}{2}},(\mathrm{100b})$

with ${\mathbf{L}}_v=\left[L_1^{-1},\dots ,L_M^{-1}\right]$ and ${\mathbf{L}}_t=\left[C_1^{-1},\dots ,C_M^{-1}\right]$.

6.4.1 Vectorized scheme for the horn equation

We write the scheme in Equations 97a, 97b in vector form by defining two $N\times 1$ length vectors $\mathit{\psi }={\left[{\psi }_0\dots {\psi }_{N-1}\right]}^T$ and $\mathit{\upsilon }={\left[{\upsilon }_{\frac{1}{2}}\dots {\upsilon }_{N+\frac{1}{2}}\right]}^T$ with spatial indices $l\in \underline{d}=0,\dots ,N-1$. The vectorized scheme is then given as

${\delta }_{t-}{\mathit{\upsilon }}^{n+\frac{1}{2}}+{\mathbf{D}}_{z+}\left({\mathbf{L}}_{\psi }{\mathit{\psi }}^n\right)+{\mathbf{G}}_p{p}_L^n+{\mathbf{C}}_v^M{\mathbf{\Phi }}^{n-1}+{\mathbf{D}}_v^M{\mu }_{t-}\left({\mathbf{L}}_{\upsilon }{\mathit{\upsilon }}^{n+\frac{1}{2}}\right)=0,(\mathrm{101a})$

${\delta }_{t+}{\mathit{\psi }}^n+{\mathbf{D}}_{z-}\left({\mathbf{L}}_{\upsilon }^{n+\frac{1}{2}}\right)\mathit{\upsilon }+{\mathbf{G}}_u{u}_{-\frac{1}{2}}^{n+\frac{1}{2}}+{\mathbf{C}}_t^M{\mathit{Q}}^{n+\frac{1}{2}}+{\mathbf{D}}_t^M{\mu }_{t+}\left({\mathbf{L}}_{\psi }{\mathit{\psi }}^n\right)=0.(\mathrm{101b})$

${\mathbf{D}}_{z+}=-{\mathbf{D}}_{z-}^T$ are the matrix forms of the spatial difference operators (Bilbao, 2009). ${\mathbf{L}}_{\psi }$ and ${\mathbf{L}}_{\upsilon }$ are $N\times N$ diagonal matrices associated with the quadratic energy densities along the bore

${\mathbf{L}}_{\psi }=\text{diag}\left(\left[\frac{\rho c^2}{S_0},\dots ,\frac{\rho c^2}{S_N-1}\right]\right),{\mathbf{L}}_{\upsilon }=\text{diag}\left(\left[\frac{S_{\frac{1}{2}}}{\rho },\dots ,\frac{S_{N+\frac{1}{2}}}{\rho }\right]\right).(\mathrm{102a})$

${\mathbf{G}}_p$ and ${\mathbf{G}}_u$ are $N\times 1$ vectors which incorporate the boundary elements by completing the spatial discretization scheme and are defined as:

${\mathbf{G}}_p={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill \cdots \hfill & \hfill 1\hfill \end{array}\right]}^T,{\mathbf{G}}_u={\left[\begin{array}{ccc}\hfill -1\hfill & \hfill \cdots \hfill & \hfill 0\hfill \end{array}\right]}^T.(\mathrm{102b})$

Viscous and thermal losses are included by means of the explicit update matricies of Equations 99a, 99b, Equations 100a, 100b. The concatenated $MN\times 1$ length viscothermal loss state vectors are $\mathit{Q}={\left[{\mathit{q}}_{\frac{1}{2}}\dots {\mathit{q}}_{N+\frac{1}{2}}\right]}^T$ and $\mathbf{\Phi }={\left[{\mathit{\varphi }}_0\dots {\mathit{\varphi }}_{N-1}\right]}^T$. The associated viscothermal loss matrices ${\mathbf{C}}_{v,t}$ and ${\mathbf{D}}_{v,t}$ are $MN\times MN$ and $N\times N$ block diagonal matrices, respectively

${\mathbf{C}}_v^M=\text{diag}\left(\left[{\mathbf{C}}_{v,\frac{1}{2}},\dots ,{\mathbf{C}}_{v,N+\frac{1}{2}}\right]\right),{\mathbf{D}}_v^M=\text{diag}\left(\left[{\mathbf{D}}_{v,\frac{1}{2}},\dots ,{\mathbf{D}}_{v,N+\frac{1}{2}}\right]\right),(\mathrm{102c})$

${\mathbf{C}}_t^M=\text{diag}\left(\left[{\mathbf{C}}_{t,0},\dots ,{\mathbf{C}}_{v,N-1}\right]\right),{\mathbf{D}}_t^M=\text{diag}\left(\left[{\mathbf{D}}_{t,0},\dots ,{\mathbf{D}}_{t,N-1}\right]\right),(\mathrm{102d})$

where the matrices at each spatial location correspond to the matrices ${\mathbf{C}}_d$ and ${\mathbf{D}}_d$ from the discrete gradient discretization procedure. An explicit form is derived by solving the linear equations in Equations 101a, 101b for ${\mathit{\upsilon }}^{n+\frac{1}{2}}$ and. ${\mathit{\psi }}^{n+1}$

${\mathit{\upsilon }}^{n+\frac{1}{2}}={\mathbf{H}}_u^{-1}\left({\mathbf{A}}_u{\mathit{\upsilon }}^{n-\frac{1}{2}}-{\mathrm{\Delta }}_t{\mathbf{D}}_{z+}{\mathbf{L}}_{\psi }{\mathit{\psi }}^n-{\mathrm{\Delta }}_t{\mathbf{C}}_v^M{\mathbf{\Phi }}^{n-1}+{\mathrm{\Delta }}_t{\mathbf{G}}_p{p}_L^n\right),(\mathrm{103a})$

${\mathit{\psi }}^{n+1}={\mathbf{H}}_p^{-1}\left({\mathbf{A}}_p{\mathit{\psi }}^n-{\mathrm{\Delta }}_t{\mathbf{D}}_{z-}{\mathbf{L}}_{\upsilon }{\mathit{\upsilon }}^{n+\frac{1}{2}}-{\mathrm{\Delta }}_t{\mathbf{C}}_t^M{\mathit{Q}}^{n+\frac{1}{2}}+{\mathrm{\Delta }}_t{\mathbf{G}}_u{u}_{-\frac{1}{2}}^{n+\frac{1}{2}}\right),(\mathrm{103b})$

with

${\mathbf{H}}_u=\mathbf{I}+\frac{1}{2}{\mathbf{D}}_v{\mathbf{L}}_{\upsilon },{\mathbf{H}}_p=\mathbf{I}+\frac{1}{2}{\mathbf{D}}_t{\mathbf{L}}_{\upsilon },(\mathrm{103c})$

${\mathbf{A}}_u=\mathbf{I}-\frac{1}{2}{\mathbf{D}}_v{\mathbf{L}}_{\psi },{\mathbf{A}}_p=\mathbf{I}-\frac{1}{2}{\mathbf{D}}_t{\mathbf{L}}_{\psi }.(\mathrm{103d})$

After computing Equations 103a, 103b, the internal viscothermal loss states ${\mathbf{\Phi }}^n$ and ${\mathit{Q}}^n$ are updated via

${\mathbf{\Phi }}^n={\mathbf{A}}_v^M{\mathbf{\Phi }}^{n-1}+\frac{1}{2}{\mathbf{B}}_v^M{\mathbf{L}}_{\upsilon }\left({\mathit{\upsilon }}^{n+\frac{1}{2}}+{\mathit{\upsilon }}^{n-\frac{1}{2}}\right),(\mathrm{104a})$

${\mathit{Q}}^{n+\frac{1}{2}}={\mathbf{A}}_t^M{\mathit{Q}}^{n+\frac{1}{2}}+\frac{1}{2}{\mathbf{B}}_t^M{\mathbf{L}}_{\psi }\left({\mathit{\psi }}^{n+1}+{\mathit{\psi }}^n\right),(\mathrm{104b})$

where—similar to ${\mathbf{C}}_{v,t}^M$ and ${\mathbf{D}}_{v,t}^M$—${\mathbf{A}}_{v,t}^M$ and ${\mathbf{B}}_{v,t}^M$ are block diagonal matrices but with dimensions $MN\times MN$ and $MN\times N$, respectively. Each diagonal entry corresponds to the update matrices at the corresponding spatial index.

The scheme in Equation 103b is driven by an external volume velocity at the input end, $l=0$ and pressure at the output end, $l=N$. $u_{-\frac{1}{2}}$ is a virtual grid point and can be determined via the equivalence ${\mu }_{t+}{u}_{-\frac{1}{2}}=u_0$. Thus,

$u_{-\frac{1}{2}}=2u_0-\frac{S_{\frac{1}{2}}}{\rho }{\upsilon }_{\frac{1}{2}}.(105)$

The pressure at the output end depends on the boundary condition at the end of the bore. This is discussed in Section 6.5 where we couple the bore to an unflanged radiation condition and in Section 7 where we couple segments of bores together through two-port elements.

Assuming that the matrix inverses are computed and applied beforehand, the lossy horn equation in Equations 103a, 103b, 103c, 103d has a computational cost of four $(N\times N)\times (N\times 1)$ matrix-vector products, two $(N\times MN)\times (MN\times 1)$ matrix-vector products, two $N\times 1$ vector-scalar products, and eight $N\times 1$ vector additions. Updating the internal loss states in Equations 104a, 104b further incurs two $(MN\times M\cdot N)\times (MN\times 1)$ matrix-vector products, two $(MN\times N)\times (N\times 1)$ matrix-vector products, two $N\times 1$ vector additions, and two $M\cdot N\times 1$ vector additions.

6.4.2 Discrete power balance

Aside from the characterization of losses as a PHS, the scheme described here is exactly the same as the scheme in Bilbao and Harrison (2016). The PHS loss subsystems can be rewritten in a manner that is symbolically equivalent to the losses in Bilbao and Harrison (2016) even though, numerically, the scheme proposed here and the scheme in Bilbao and Harrison (2016) are not equivalent. The discrete energetic quantities follow from the proof given in that article. As such, we do not repeat the proof here. The horn equation scheme is conditionally stable and is determined by the Courant-Friedrichs-Lewy (CFL) condition $\lambda \le 1$ with $\lambda =\frac{c{\mathrm{\Delta }}_t}{{\mathrm{\Delta }}_z}$ (Bilbao and Harrison, 2016). The discrete Hamiltonian associated with Equation 97 is,

${\mathfrak{h}}_b={\mathfrak{h}}_p+{\mathfrak{h}}_{u,mod}+{\mathfrak{h}}_t+{\mathfrak{h}}_{v,mod}(106)$

with

${\mathfrak{h}}_p=\frac{1}{2}{\left({‖\sqrt{L_{\psi }}\psi ‖}_d^{\prime }\right)}^2,\hspace{1em}{\mathfrak{h}}_{u,mod}=\frac{1}{2}{⟨\sqrt{L_{\upsilon }}\upsilon ,e_{t-}\upsilon ⟩}_{\underline{d}},(\mathrm{107a})$

${\mathfrak{h}}_t=\frac{1}{2}{\left({‖{\mathit{q}}^T{\mathbf{L}}_t\mathit{q}‖}_d^{\prime }\right)}^2,\hspace{1em}{\mathfrak{h}}_{v,mod}=\frac{1}{2}{⟨{\mathbf{L}}_v\mathit{\varphi },e_{t-}\mathit{\varphi }⟩}_{\underline{d}},(\mathrm{107b})$

where $\Vert \cdot {\Vert }_{\underline{d}}$ and ${⟨\cdot ,\cdot ⟩}_{\underline{d}}$ denote the norm and inner product over the spatial domain $\underline{d}=\mathrm{0,1},\dots ,N-1$. Conversely, $\Vert \cdot {\Vert }_d^{\prime }$ and ${⟨\cdot ,\cdot ⟩}_d^{\prime }$ denote the primed norm and inner product over the entire domain $d=\mathrm{0,1},\dots ,N$. Readers are referred to Bilbao (2009) regarding the definition of these operators. The dissipated energies are

${\mathfrak{q}}_t={⟨{\delta }_{t+}\mathit{q},{\mathbf{R}}_t{\delta }_{t+}\mathit{q}⟩}_d^{\prime },(\mathrm{108a})$

${\mathfrak{q}}_{v,mod}={⟨\widehat{\mathit{v}},{\mathbf{R}}_v{\mu }_{t+}{\mu }_{t-}\widehat{\mathit{v}}⟩}_{\underline{d}},{\mathfrak{q}}_{v_0,mod}={⟨L_{\upsilon }\upsilon ,R_0{\mu }_{t+}{\mu }_{t-}{L}_{\upsilon }\upsilon ⟩}_{\underline{d}},(\mathrm{108b})$

with

$\widehat{\mathit{v}}=\mathbf{G}\upsilon -{\mathbf{L}}_v\mathit{\varphi },(\mathrm{108c})$

and the discrete power transferred through the boundaries is,

${\mathfrak{b}}_b=\left({\mu }_{t+}{p}_N^n\right)u_N^{n+\frac{1}{2}}-\left({\mu }_{t+}{p}_0^n\right)u_0^{n+\frac{1}{2}}.(109)$

Due to discretization with the Störmer-Verlet method, the discrete stored energy and dissipated energy is not an exact sampling of the continuous energy function. Quantities denoted with the subscript $mod$ are only guaranteed non-negative if $\lambda \le 1$ and the system instead maintains a modified energy balance. Which energies are modified depends on the alignment with the boundary variables, in this case the modified energies are related to $v$ and $p_v$ are centered on time $n$ and the interleaved spatial grid $l+\frac{1}{2}$.

6.5 Radiation model

We use the second-order circuit model, shown in Figure 6, presented by Bilbao and Chick (2013) to model the radiation from an unflanged end of the bore. The model is described by the following PHS

$\underset{{\mathbf{f}}_{\mathit{rad}}}{\underbrace{\left[\begin{array}{c}\hfill {\partial }_t{\varphi }_{\mathit{rad}}\hfill \\ \hfill {\partial }_t{q}_{\mathit{rad}}\hfill \end{array}\right]}}=\left((\right.\left[\begin{array}{cc}\hfill 0\hfill & \hfill -1\hfill \\ \hfill 1\hfill & \hfill 0\hfill \end{array}\right]-\underset{{\mathbf{R}}_{\mathit{rad}}}{\underbrace{\left[\begin{array}{cc}\hfill R_{\mathit{rad,1}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill R_{\mathit{rad,2}}^{-1}\hfill \end{array}\right]}}\left.\right))\underset{{\mathbf{e}}_{\mathit{rad}}}{\underbrace{\left[\begin{array}{c}\hfill {\partial }_{{\varphi }_{\mathit{rad}}}{\mathcal{H}}_{\mathit{rad}}\hfill \\ \hfill {\partial }_{q_{\mathit{rad}}}{\mathcal{H}}_{\mathit{rad}}\hfill \end{array}\right]}}+\left((\right.\underset{{\mathbf{G}}_{\mathit{rad}}}{\underbrace{\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 1\hfill \end{array}\right]}}-\underset{{\mathbf{P}}_{\mathit{rad}}}{\underbrace{\left[\begin{array}{c}\hfill R_{\mathit{rad,1}}\hfill \\ \hfill 0\hfill \end{array}\right]}}\left.\right))u_L,(\mathrm{110a})$

$p_L={\left({\mathbf{G}}_{\mathit{rad}}+{\mathbf{P}}_{\mathit{rad}}\right)}^T{\mathbf{e}}_{\mathit{rad}}+R_{\mathit{rad,1}}{u}_L,(\mathrm{110b})$

with parameters

$R_{\mathit{rad,1}}=Z_{c,L},R_{\mathit{rad,2}}=0.505Z_{c,L},L_{\mathit{rad}}=0.613\frac{r}{c}{Z}_{c,L},C_{\mathit{rad}}=1.111\frac{r}{c}\frac{1}{Z_{c,L}},(\mathrm{110c})$

where $Z_{c,L}=\frac{\rho c}{S_L}$ the characteristic impedance at the end of the bore. The Hamiltonian of the unflanged radiation system is:

${\mathcal{H}}_{\mathit{rad}}=\frac{{\varphi }_{\mathit{rad}}^2}{2L_{\mathit{rad}}}+\frac{q_{\mathit{rad}}^2}{2C_{\mathit{rad}}},(111)$

and the dissipated energy and boundary power are

${\mathcal{Q}}_{\mathit{rad}}=⟨{\left[{\mathbf{e}}_{\mathit{rad}}{u}_L\right]}^T,\left[\begin{array}{cc}\hfill {\mathbf{R}}_{\mathit{rad}}\hfill & \hfill {\mathbf{P}}_{\mathit{rad}}\hfill \\ \hfill {\mathbf{P}}_{\mathit{rad}}^T\hfill & \hfill R_{\mathit{rad,1}}\hfill \end{array}\right]{\left[{\mathbf{e}}_{\mathit{rad}}{u}_L\right]}^T⟩,{\mathcal{B}}_{\mathit{rad}}=p_L{u}_L.(112)$

Figure 6

Figure 6. Circuit approximation of the radiation from an unflanged end of a circular bore after Bilbao and Chick (2013)

We discretize this system with the discrete gradient method on the $n+\frac{1}{2}$ time grid and produce the explicit state update and output equations

${\mathbf{x}}_{\mathit{rad}}^{n+1\frac{1}{2}}={\mathbf{A}}_r{\mathbf{x}}_{\mathit{rad}}^{n+\frac{1}{2}}+{\mathbf{B}}_r{u}_L^{n+\frac{1}{2}},(\mathrm{113a})$

$p_L^{n+\frac{1}{2}}={\mathbf{C}}_r{\mathbf{x}}_{\mathit{rad}}^{n+\frac{1}{2}}+D_r{u}_L^{n+\frac{1}{2}}.(\mathrm{113b})$

The discrete energy balance is maintained by the exact discretized forms of the quantities in Equations 111, 112. We will discuss coupling the radiation to the bore in Section 8.

7 Tonehole model

The woodwind tonehole is normally represented by a frequency-domain equivalent circuit model involving four elements: two series inertances equal to $Z_a/2$, an inner length correction impedance $Z_i$, and a shunt impedance $Z_s$ as shown on the left of Figure 7 (Chaigne and Kergomard, 2016; Keefe, 1982; Dubos et al., 1999; Dalmont et al., 2002; Lefebvre and Scavone, 2012). The values corresponding to the lumped circuit elements are usually obtained via a fit to 3D simulation data of geometries with sideholes (Dubos et al., 1999; Lefebvre and Scavone, 2012) or analytical Green’s function approaches (Keefe, 1982; Dubos et al., 1999). The one-dimensional model has also been validated against measurements in Dalmont et al. (2002). In this section, we will first present the lumped model parameters based on (Lefebvre and Scavone, 2012). We will then propose a low-frequency circuit approximations of the open and closed tonehole shunt immitances before uniting the models through a switching PHS. The switching PHS model is then discretized with the discrete gradient method.

Figure 7

Figure 7. Left: Equivalent circuit model of the tonehole with negative series inertances $Z_a/2$, inner length correction impedance $Z_i$, and shunt impedance $Z_s$. Right: Low-frequency switching model used in this article with switch $s$. $Z_{\mathit{rth}}$ corresponds to the radiation impedance model given in Section 6.5 dependent on the tonehole radius.

7.1 Lumped tonehole parameters

For a tonehole with tonehole radius $b$, tonehole height $t_h$, and main bore radius $r$, the series inertance $Z_a$ is described by

$Z_a=j\omega \frac{\rho }{S}{t}_a^{\left(o,c\right)},(\mathrm{114a})$

which has different values when the tonehole is open and closed corresponding to the length correction values $t_a^{(o)}$ and $t_a^{(c)}$:

$t_a^{\left(o\right)}=-b{d}^2\left(0.36-0.06\tanh \left(2.7t_h/b\right)\right),t_a^{\left(c\right)}=-b{d}^2\left(0.12+0.17\tanh \left(2.4t_h/b\right)\right),(\mathrm{114b})$

with $d=b/r$. For the open tonehole, the lumped model shunt impedance is equal to:

$Z_s^{\left(o\right)}=j\omega \frac{\rho c}{S_h}{t}_i+j\frac{\rho c}{S_h}\tan \left[\frac{\omega }{c}\left(t_h+t_m+t_r\right)\right].(115)$

$S_h=\pi b^2$ is the cross-sectional area of the tonehole. $t_i$, $t_m$, and $t_r$ are length correction parameters associated with the inner length correction, the matching volume correction, and the radiation at the end of the tonehole. The matching volume and radiation length correction are

$t_m=\frac{bd}{8}\left(1+0.207d^3\right),t_r=\arctan \left(\frac{Z_{\mathit{rth}}}{j\rho c/S_h}\right),(116)$

and the inner length correction is

$t_i=\left(0.822-0.095d-1.566d^2+2.138d^3-1.640d^4+0.502d^5\right)b.(\mathrm{117a})$

$Z_{\mathit{rth}}$ is the radiation impedance associated with the end condition of the open tonehole. Lefebvre and Scavone (2012) propose a frequency dependent length correction such that $t_i(k)=t_i\times G(d,kb)$ with $kb=\omega b/c$,

$G\left(d,kb\right)=1+H\left(d\right)I\left(kb\right),(\mathrm{117b})$

$H\left(d\right)=1-4.56d+6.55d^2,I\left(kb\right)=0.17ka+0.92{\left(kb\right)}^2+0.16{\left(kb\right)}^3-0.29{\left(kb\right)}^4.(\mathrm{117c})$

The closed tonehole impedance includes the same matching volume and inner length correction terms as appear in the open tonehole impedance. The form of the closed tonehole impedance is different due to the closed end condition and, correspondingly, the radiation length correction is omitted,

$Z_s^{\left(c\right)}=j\omega \frac{\rho c}{S_h}{t}_i-j\frac{\rho c}{S_h}\cot \left[\frac{\omega }{c}\left(t_h+t_m\right)\right].(118)$

In both models, viscothermal losses can be incorporated by replacing the imaginary wavenumber $jk=j\omega /c$ with the lossy wavenumber $\mathrm{\Gamma }$ (Dalmont et al., 2002). $\mathrm{\Gamma }$ is related to the viscothermal losses by:

$\mathrm{\Gamma }={\left[\left(j\omega \frac{\rho }{S_h}+Z_v\left(b\right)\right)\left(j\omega \frac{S_h}{\rho c^2}+Y_t\left(b\right)\right)\right]}^{1/2},(119)$

where the related loss immitances are evaluated with the tonehole radius $b$.

7.2 Low-frequency tonehole model

A known hindrance in developing a time-domain model based on the lumped tonehole approximation is that the series length correction, $t_a$, is always negative whether the tonehole is open or closed. This represents a slowing of the flow as the bore interacts with the tonehole (Keefe, 1982). It is likely that the contribution of the negative inertance should not be ignored. We take the approach in van Walstijn and Scavone (2000) and subsume the negative length into the length of the bore. Including the negative length in this manner has the effect of “breaking” the woodwind bore into a system of segments connected by two-port junctions. We will not model the effect of changing the bore length to account for $t_a$ when the tonehole is open and closed. Instead, we will use the open tonehole $t_a^{(o)}$ in all cases. This approximation is motivated by the fact that the difference in these length corrections is minimal and not the dominant effect of the tonehole. The resulting low-frequency approximations of the open shunt impedance and closed shunt admittance compared to the lossless and lossy tonehole models from (Lefebvre and Scavone, 2012) is displayed in Figure 8. The responses are plotted against the related frequency values of $kb=\omega /cb$ up to the cutoff frequency $(kb=1.84)$ of the planar tonehole mode. Above this cutoff frequency, the tonehole model is no longer valid (Chaigne and Kergomard, 2016, 7.6.2).

Figure 8

Figure 8. Comparison of the normalized open tonehole shunt impedance (A) and closed tonehole shunt admittance (B) for the lossless (red, dashed), lossy (blue, solid) tonehole models in Lefebvre and Scavone (2012), and the circuit approximation (green, solid) proposed in this article. Parameters used for the plots are $r=7.5$mm, $\delta =0.5$, $t_h=1.1b$, $T$ = 26.85$°$C, $kb=\omega /cb$.

7.2.1 Open shunt impedance

Following the approach in van Walstijn and Scavone (2000), we develop a discrete-time model of the closed and open shunt impedances via a low-frequency approximation. To model the open shunt impedance, we use the frequency independent $t_i$ in Equation 117a and a small angle approximation for $\tan (\cdot )$ in Equation 115. We cannot use the frequency dependent form of $t_i$ because the function is not positive real and cannot be represented by any circuit network. Radiation at the end of the tonehole is modeled using the circuit representation presented in Section 6.5. The low-frequency open shunt impedance is then given by the series combination of two inertances and the impedance of the unflanged radiation approximation

$Z_{s,lf}^{\left(o\right)}=j\omega L_i+j\omega L_o+Z_{\mathit{rth}}\left(b\right),(\mathrm{120a})$

$L_i=\frac{\rho }{S_h}{t}_i,L_o=\frac{\rho }{S_h}\left(t_h+t_m\right).(\mathrm{120b})$

$Z_{\mathit{rth}}$ depends on the tonehole radius, not the bore radius. In Figure 8, the behavior of the circuit approximation given by Equations 120a, 120b is compared to the normalized open tonehole shunt impedance in Lefebvre and Scavone (2012) with and without losses. The low-frequency approximation is a good match to both the lossy and lossless theoretical response for $kb<0.5$.

7.2.2 Closed shunt admittance

For the closed tonehole, the small angle approximation of $\cot (\cdot )$ results in a compliance in series with an inductance. In our model, we include a fictitious resistance $R_c$ in series with our compliance. The resulting closed tonehole model is a series RLC circuit as was presented in Section 4. The additional resistance is included to better match the effects of viscothermal losses and the frequency dependent inner length correction in the closed tonehole model.

In the absence of these components, the closed tonehole model in Equation 118 behaves akin to a Helmholtz resonator and will produce a strong resonance at a frequency proportional to the volume of the closed tonehole (Chaigne and Kergomard, 2016, 1.5). One effect of the viscothermal losses (Equation 119) and the frequency dependent inner length correction (Equations 117a, 117b, 117c) is to dramatically dampen the tonehole resonance. A secondary effect is a shifting down—in frequency—of the same resonance. In both cases, the frequency dependent length correction is the more dominant effect. The additional resistance, $R_c$, allows us to tune the strength of the resonance in our low-frequency approximation to better match the peak admittance in the theoretical model. The compliance, $C_c$ can also be adjusted to account for the shift in tonehole resonance. This is shown in Figure 8 where the both parameters, $R_c$ and $C_c$, have been adjusted to fit the model from Lefebvre and Scavone (2012).

Both parameters can also be used to compensate for the effects of frequency warping. When the tonehole model is eventually discretized, the theoretical closed tonehole resonance can potentially be at a frequency near or above the Nyquist limit. To minimize the effects of frequency warping, $R_c$ can be adjusted so that the peak admittance in the low-frequency model is equal to the peak admittance below the Nyquist limit. $C_c$ can also be adjusted to prewarp the response.

7.3 Switching PHS tonehole model

To transition between the closed and open tonehole models, van Walstijn and Scavone (2000) proposed placing the shunt compliance and inertance, corresponding to the closed and open tonehole models, in parallel. The compliance and inertances values are then modified as a function of a variable between zero and one. The overall effect is to open the circuit on the opposing parallel branch by setting the compliance to zero when the tonehole is open or the inertance to infinity when the tonehole is closed. An immediate issue with this approach, outside of simulating a system with infinite inertance at times, is that the energy stored in the tonehole system becomes infinite when the tonehole is opened.

In light of these issues, we instead model the tonehole transition using a switching PHS which is used to represent systems with variable topology (Duindam et al., 2009, 2.2.5). As a result, the energy stored in the tonehole is independent of the tonehole state. The switch does not have to be instantaneous and half-holing is still possible using the switching model. We combine the open and closed tonehole systems using the switching variable $s$. When $s=0$, the tonehole is completely closed and when $s=1$ the tonehole is completely open. The PHS formulation of the switching tonehole system is given below, and the complete model is represented as a circuit in Figure 7

$\left[\begin{array}{c}\hfill v_i={\partial }_t{\varphi }_i\hfill \\ \hfill i_c={\partial }_t{q}_c\hfill \\ \hfill i_o={\partial }_{{\varphi }_o}{\mathcal{H}}_{th}\hfill \\ \hfill v_{\mathit{rth}}={\partial }_t{\varphi }_{\mathit{rth}}\hfill \\ \hfill i_{\mathit{rth}}={\partial }_t{q}_{r,th}\hfill \end{array}\right]=\left({\mathbf{J}}_{th}\left(s\right)-{\mathbf{R}}_{th}\left(s\right)\right)\underset{{\mathbf{e}}_{th}}{\underbrace{\left[\begin{array}{c}\hfill i_i={\partial }_{{\varphi }_i}{\mathcal{H}}_{th}\hfill \\ \hfill v_c={\partial }_{q_c}{\mathcal{H}}_{th}\hfill \\ \hfill v_o={\partial }_t{\varphi }_o\hfill \\ \hfill i_{r,th}={\partial }_{{\varphi }_{\mathit{rth}}}{\mathcal{H}}_{th}\hfill \\ \hfill v_{r,th}={\partial }_{q_{\mathit{rth}}}{\mathcal{H}}_{th}\hfill \end{array}\right]}}+{\mathbf{G}}_{th}\left[\begin{array}{c}\hfill {\stackrel{-}{u}}_L\hfill \\ \hfill {\stackrel{+}{p}}_0\hfill \end{array}\right],(\mathrm{121a})$

$\left[\begin{array}{c}\hfill {\stackrel{-}{p}}_L\hfill \\ \hfill {\stackrel{+}{u}}_0\hfill \end{array}\right]={\left({\mathbf{G}}_{th}\right)}^T{\mathbf{e}}_{th}+\underset{{\mathbf{M}}_{th}}{\underbrace{\left[\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 0\hfill \end{array}\right]}}\left[\begin{array}{c}\hfill {\stackrel{-}{u}}_L\hfill \\ \hfill {\stackrel{+}{p}}_0\hfill \end{array}\right],(\mathrm{121b})$

with

${\mathbf{J}}_{th}\left(s\right)=\left[\begin{array}{ccccc}\hfill 0\hfill & \hfill s-1\hfill & \hfill -s\hfill & \hfill 0\hfill & \hfill -s\hfill \\ \hfill 1-s\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill s\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -s\hfill \\ \hfill s\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill s\hfill & \hfill 0\hfill \end{array}\right],(\mathrm{122a})$

${\mathbf{R}}_{th}\left(s\right)=\left[\begin{array}{ccccc}\hfill s{R}_{\mathit{rth,1}}+\left(1-s\right)R_c\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill s{R}_{\mathit{rth,1}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill s{R}_{\mathit{rth,1}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill s{R}_{\mathit{rth,1}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill s{R}_{\mathit{rth,2}}^{-1}\hfill \end{array}\right],(\mathrm{122b})$

${\mathbf{G}}_{th}={\left[\begin{array}{ccccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T.(\mathrm{122c})$

The overset symbols $\stackrel{-}{()}$ and $\stackrel{+}{()}$ denote the acoustic variables to the left and right of the tonehole, respectively. The total stored energy in the tonehole is

${\mathcal{H}}_{th}=\frac{{\varphi }_i^2}{2L_i}+\frac{q_c^2}{2C_c}+\frac{{\varphi }_o^2}{2L_o}+\frac{{\varphi }_{\mathit{rth}}^2}{2L_{\mathit{rth}}}+\frac{q_{\mathit{rth}}^2}{2C_{\mathit{rth}}}.(123)$

The energy dissipated in the system is now also a function of $s$,

${\mathcal{Q}}_{th}\left(s\right)=⟨{\mathbf{e}}_{th},{\mathbf{R}}_{th}\left(s\right){\mathbf{e}}_{th}⟩,(124)$

and the power transmitted through the two-port boundary is

${\mathcal{B}}_{th}={\stackrel{-}{p}}_L{\stackrel{-}{u}}_L+{\stackrel{+}{p}}_0{\stackrel{+}{u}}_0.(125)$

The model is a PHS for the tonehole state value $0\le s\le 1$ as evident by the skew-symmetric structure of ${\mathbf{J}}_{th}(s)$. In Equations 121a, 121b, we have a hybrid definition of efforts and flows due to the open shunt inertance, $L_o$, and the inner length correction inertance, $L_i$, sharing an effort. Typically, these two inertances would be combined into a single element to eliminate the non-causal connection (Hélie, 2022). However, the two elements must remain distinct so that $L_i$ can interact with the elements on the closed tonehole branch of the model. A similar approach is taken in (Müller and Hélie, 2021; Section 6.2) to keep inductances distinct.

7.3.1 Discretization of the switching tonehole model

The system in Equations 121a, 121b has a hybrid PHS definition that is not easily discretized using the methods outlined in Section 3.2. However, using the definition in Equation 13 we can rewrite Equation 121a as

${\delta }_{t+}\left({\mathbf{L}}_{f_1}{\mathbf{x}}_{th}^{n+\frac{1}{2}}\right)+{\mathbf{L}}_{f_2}{\mathbf{x}}_{th}^{n+\frac{1}{2}}=\left({\mathbf{J}}_{th}\left(s\right)-{\mathbf{R}}_{th}\left(s\right)\right)\left({\mu }_{t+}\left({\mathbf{L}}_{e_1}{\mathbf{x}}_{th}^{n+\frac{1}{2}}\right)-{\mathbf{L}}_{e_2}{\mathbf{x}}_{th}^{n+\frac{1}{2}}\right)+{\mathbf{G}}_{th}{\mathbf{u}}_{th}^{n+\frac{1}{2}},(126)$

where ${\mathbf{x}}_{th}={\left[\begin{array}{ccccc}\hfill {\varphi }_i\hfill & \hfill q_c\hfill & \hfill {\varphi }_o\hfill & \hfill {\varphi }_{\mathit{rth}}\hfill & \hfill q_{\mathit{rth}}\hfill \end{array}\right]}^T$ are the tonehole energy storage states and the system is discretized on the $n+\frac{1}{2}$ timestep. The matrices denoted with $\mathbf{L}$ are diagonal matrices equal to

${\mathbf{L}}_{f_1}=\text{diag}\left(\left[\mathrm{1,1},{\mathrm{\Delta }}_t/\left(2L_o\right),\mathrm{1,1}\right]\right),(\mathrm{127a})$

${\mathbf{L}}_{f_2}=\text{diag}\left(\left[\mathrm{0,0},{\mathrm{\Delta }}_t/\left(2L_o\right),\mathrm{0,0}\right]\right),(\mathrm{127b})$

${\mathbf{L}}_{e_1}=\text{diag}\left(\left[1/L_i,1/C_c,2/{\mathrm{\Delta }}_t,1/L_{\mathit{rth}},1/C_{\mathit{rth}}\right]\right),(\mathrm{128a})$

${\mathbf{L}}_{e_2}=\text{diag}\left(\left[\mathrm{0,0,2}/{\mathrm{\Delta }}_t,\mathrm{0,0}\right]\right).(\mathrm{128b})$

The explicit update for the tonehole system is

${\mathbf{x}}_{th}^{n+1\frac{1}{2}}={\mathbf{A}}_{th}{\mathbf{x}}_{th}^{n+\frac{1}{2}}+{\mathbf{B}}_{th}\left[\begin{array}{c}\hfill {\stackrel{-}{u}}_L^{n+\frac{1}{2}}\hfill \\ \hfill {\stackrel{+}{p}}_0^{n+\frac{1}{2}}\hfill \end{array}\right],(\mathrm{129a})$

$\left[\begin{array}{c}\hfill {\stackrel{-}{p}}_L^{n+\frac{1}{2}}\hfill \\ \hfill {\stackrel{+}{u}}_0^{n+\frac{1}{2}}\hfill \end{array}\right]={\mathbf{C}}_{th}{\mathbf{x}}_{th}^{n+\frac{1}{2}}+{\mathbf{D}}_{th}\left[\begin{array}{c}\hfill {\stackrel{-}{u}}_L^{n+\frac{1}{2}}\hfill \\ \hfill {\stackrel{+}{p}}_0^{n+\frac{1}{2}}\hfill \end{array}\right],(\mathrm{129b})$

and defined by the matrices

${\mathbf{H}}_{th}={\left({\mathbf{L}}_{f_1}-\frac{{\mathrm{\Delta }}_t}{2}\mathbf{A}{\mathbf{L}}_{e_1}\right)}^{-1},{\mathbf{A}}_{th}={\mathbf{H}}_{th}\left[{\mathbf{L}}_{f_1}+{\mathrm{\Delta }}_t\left(\frac{1}{2}\mathbf{A}{\mathbf{L}}_{e_1}-\mathbf{A}{\mathbf{L}}_{e_2}-{\mathbf{L}}_{f_2}\right)\right],(\mathrm{129c})$

${\mathbf{B}}_{th}={\mathrm{\Delta }}_t{\mathbf{H}}_{th}{\mathbf{G}}_{th},(\mathrm{129d})$

${\mathbf{C}}_{th}={\mathbf{G}}_{th}^T\left(\frac{1}{2}{\mathbf{L}}_{e_1}\left({\mathbf{A}}_{th}+\mathbf{I}\right)+{\mathbf{L}}_{e_2}\right),{\mathbf{D}}_{th}=\frac{1}{2}{\mathbf{G}}_{th}^T{\mathbf{L}}_{e_1}{\mathbf{B}}_{th}+{\mathbf{M}}_{th},(129e)$

where $\mathbf{A}=\left({\mathbf{J}}_{th}(s)-{\mathbf{R}}_{th}(s)\right)$ and $\mathbf{I}$ is the appropriately sized $5\times 5$ identity matrix. All the matrices are dependent on the tonehole state, $s$, and must be recomputed if the tonehole state changes. The discrete energy is an exact discretization of the terms in Equations 123–125.

In Figure 9, we examine the input impedance of a cylindrical pipe with a single tonehole at its center for different values of $s$. We model the system using the discrete PHS described thus far. The simulation includes losses and the unflanged radiation condition at the end of the bore. In the case the tonehole is completely closed $(s=0)$ or open $(s=1)$, the response is compared to the input impedance derived using the transfer matrix method (TMM) (Scavone, 2024). The discrete PHS model closely matches the behavior of the TMM at frequencies up to around 16 kHz. As $s$ goes from 0 to 1 the response gradually shifts, starting with the higher frequencies.

Figure 9

Figure 9. Input impedances of a cylindrical tube and single tonehole derived from a discrete PHS simulation (blue) for various tonehole conditions $s$. The cylindrical tube has length $L=20f_s/c-t_a$ m and radius $r=7.5$mm with a tonehole at $L/2$ with parameters $\delta =0.5$ and $t_h=1.1b$. The simulation was carried out with $T=26.85°$C and $f_s=48$kHz. At the fully closed $(s=0)$ and open $(s=1)$ states the response is compared to the input impedance derived using the TMM (red, dashed).

8 A complete model

To build a complete model we must connect the bore, tonehole, radiation, and reed sub-models discussed thus far. The pressure and volume velocities at the external ports in these elements have been defined such that they are linked by the $2\times 2$ skew-symmetric matrix, $\mathbf{J}$, defining a parallel connection. The pressure output from one model at $z=l$ is equal to the pressure going into another model $(p_l^{\text{out}}=p_l^{\text{in}})$ and the volume velocities have opposite sign $(u_l^{\text{out}}=-u_l^{\text{in}})$. We will first discuss coupling two bores together through a general two-port system and coupling a bore to the radiation model in Section 6.5. We will then discuss the non-linear coupling between the bore and reed model from Section 5. Finally, we will perform a simulation using a simplified bore profile with toneholes.

8.1 Linear coupling to two-ports and one-port systems

We first consider a general two-port system coupling two bores discretized with the scheme in Equations 97a, 97b. We seek to solve the equation,

$\left[\begin{array}{c}\hfill {\mu }_{t+}\left({\stackrel{-}{p}}_N^n\right)\hfill \\ \hfill -{\stackrel{+}{u}}_0^{n+\frac{1}{2}}\hfill \end{array}\right]=\mathbf{C}{\mathbf{x}}^{n+\frac{1}{2}}+\mathbf{D}\left[\begin{array}{c}\hfill {\stackrel{-}{u}}_N^{n+\frac{1}{2}}\hfill \\ \hfill {\mu }_{t+}\left({\stackrel{+}{p}}_0^n\right)\hfill \end{array}\right],(130)$

for ${\stackrel{-}{p}}_N^{n+1}$ and ${\stackrel{+}{u}}_0^{n+\frac{1}{2}}$ at time $n+\frac{1}{2}$ following the first step of the Störmer-Verlet scheme. The overset symbols $\stackrel{-}{()}$ and $\stackrel{+}{()}$ will again be used denote elements pertaining to the system to the left and right of the two-port. $\mathbf{C}$ and $\mathbf{D}$ are the explicit output matrices produced by the discretization of the two-port junction. $\mathbf{x}$ refers to the internal storage states for the two-port junction.

Using Equations 98a, 98b, 13, we rewrite ${\stackrel{-}{u}}_N^{n+\frac{1}{2}}$ in terms of ${\mu }_{t+}\left({\stackrel{-}{p}}_N^n\right)$ and some additional terms. The same can be done to express ${\mu }_{t+}\left({\stackrel{+}{p}}_0^n\right)$ in terms of ${\stackrel{+}{u}}_0^{n+\frac{1}{2}}$ and some additional terms. Then we can rewrite Equation 130 as a function of a diagonal matrix, $\mathbf{V}$, and a residue vector containing the additional terms $\mathbf{r}$,

$\left[\begin{array}{c}\hfill {\mu }_{t+}\left({\stackrel{-}{p}}_N^n\right)\hfill \\ \hfill -{\stackrel{+}{u}}_0^{n+\frac{1}{2}}\hfill \end{array}\right]=\mathbf{C}{\mathbf{x}}^{n+\frac{1}{2}}+\mathbf{D}\left(\mathbf{V}\left[\begin{array}{c}\hfill {\mu }_{t+}\left({\stackrel{-}{p}}_N^n\right)\hfill \\ \hfill {\stackrel{+}{u}}_0^{n+\frac{1}{2}}\hfill \end{array}\right]+\mathbf{r}\right).(131)$

This equation can be solved for the vector ${\left[{\mu }_{t+}\left({\stackrel{-}{p}}_N^n\right){\stackrel{+}{u}}_0^{n+\frac{1}{2}}\right]}^T$ and then ${\stackrel{-}{p}}_N^{n+1}$ is found by definition of the forward averaging operator. For the lossy horn equation, the diagonal matrix $\mathbf{V}$ and the residue vector $\mathbf{r}$ are equal to

$\mathbf{V}=\text{diag}\left(\left[-\frac{{\stackrel{-}{\mathrm{\Delta }}}_z}{2{\stackrel{-}{L}}_{\psi }}\left(\frac{2}{{\mathrm{\Delta }}_t}+{\stackrel{-}{D}}_t\right),\frac{2{\stackrel{+}{L}}_{\psi }{\mathrm{\Delta }}_t}{{\stackrel{+}{\mathrm{\Delta }}}_z\left({\stackrel{+}{D}}_t{\mathrm{\Delta }}_t+2\right)}\right]\right),(\mathrm{132a})$

and

$\mathbf{r}=\left[\begin{array}{c}\hfill {\stackrel{-}{u}}_{N-\frac{1}{2}}^{n+\frac{1}{2}}+\frac{\stackrel{-}{{\mathrm{\Delta }}_z}}{2}\left[\frac{2}{{\mathrm{\Delta }}_t{\stackrel{-}{L}}_{\psi }}{\stackrel{-}{p}}_N^n-{\stackrel{-}{\mathbf{C}}}_t{\stackrel{-}{\mathit{q}}}_t^{n-\frac{1}{2}}\right]\hfill \\ \hfill \frac{{\stackrel{+}{\mathrm{\Delta }}}_z\left(2{\stackrel{+}{p}}_0^n-{\mathrm{\Delta }}_t{\stackrel{+}{L}}_{\psi }{\stackrel{+}{\mathbf{C}}}_t{\stackrel{+}{\mathit{q}}}_t^{n-\frac{1}{2}}\right)-2{\mathrm{\Delta }}_t{\stackrel{+}{L}}_{\psi }{\stackrel{+}{u}}_{\frac{1}{2}}^{n+\frac{1}{2}}}{{\stackrel{+}{\mathrm{\Delta }}}_z\left({\stackrel{+}{D}}_t{\mathrm{\Delta }}_t+2\right)}\hfill \end{array}\right].(\mathrm{132b})$

The overset symbols have again been used to distinguish elements $L_{\psi }=\frac{\rho c^2}{S_N}$ at the bore boundaries, ${\mathrm{\Delta }}_z$ the spatial discretization step for each bore, and ${\mathbf{C}}_t$, $D_t$, and ${\mathit{q}}_t^{n-\frac{1}{2}}$ corresponding to the explicit thermal loss output matrices and state at the end of each bore. After solving Equation 131, the next input values are used to update the two-port junction states and the thermal loss states at the end of the left bore. The system in Equation 131 is only a $2\times 2$ system and can be solved without the use of a matrix inverse. In the case when there is a sharp discontinuity in the bore profile, the same approach can be used to couple two bores of different radii. The discontinuity two-port is memoryless with $\mathbf{C}=\mathbf{0}$ and

$\mathbf{D}=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 0\hfill \end{array}\right].(133)$

To couple the bore to the one-port lumped radiation model described Section 6.5 we only need to solve the top line of Equation 131 and the matrix $\mathbf{D}$ reduces to a scalar.

8.2 Nonlinear coupling between the reed and bore

The final remaining step is to couple the reed model from Section 5 to the lossy horn equation at $l=0$. Without any loss of generality, we shift the discrete-time reed model to the $n+\frac{1}{2}$ time step. Then, due to the aforementioned sign conventions, the discrete-time relationship between the volume velocity at the input to the bore and the output of the reed is given as

$u_0^{n+\frac{1}{2}}=S_r{\mu }_{t+}\left(\frac{{\nu }^{n+\frac{1}{2}}}{m}\right)+\text{sign}\left(p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right)S_j^{n+\frac{1}{2}}\sqrt{\frac{2\left|p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right|}{\rho }},(134)$

corresponding to the bottom line of Equation 130. The discrete-time definitions for the jet-area $S_j$ and pressure difference $p_{\mathrm{\Delta }}$ are

$S_j^{n+\frac{1}{2}}={\left[y_l-y^{n+\frac{1}{2}}\right]}_{+},p_{\mathrm{\Delta }}^{n+\frac{1}{2}}=p_m^{n+\frac{1}{2}}-{\mu }_{t+}{p}_0^n.(135)$

The associated discrete energy loss due to the Bernoulli flow and power transfer through the boundary are

${\mathfrak{q}}_f=S_j^{n+\frac{1}{2}}\text{sign}\left(p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right)\sqrt{\frac{2}{\rho }}{\left|p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right|}^{3/2},{\mathfrak{b}}_r=p_{in}^{n+\frac{1}{2}}{u}_{in}^{n+\frac{1}{2}}+p_m^{n+\frac{1}{2}}{u}_m^{n+\frac{1}{2}}.(136)$

Using the definition of $u_0^{n+\frac{1}{2}}$ from Equation 98a and Equations 13, 135, we rewrite $u_0^{n+\frac{1}{2}}$ in terms of $p_{\mathrm{\Delta }}^{n+\frac{1}{2}}$

$u_0^{n+\frac{1}{2}}=-b_1{p}_{\mathrm{\Delta }}^{n+\frac{1}{2}}+b_1{p}_m^{n+\frac{1}{2}}+b_0^n,(137)$

where

$b_1=\frac{{\mathrm{\Delta }}_z}{2}\left(D_t+\frac{2}{L_{\psi }{\mathrm{\Delta }}_t}\right),(\mathrm{138a})$

$b_0^{n+\frac{1}{2}}=u_{\frac{1}{2}}^{n+\frac{1}{2}}+\frac{{\mathrm{\Delta }}_z}{2}\left[{\mathbf{C}}_t{\mathit{q}}^{n-\frac{1}{2}}-\frac{2}{L_{\psi }{\mathrm{\Delta }}_t}{p}_0^n\right].(\mathrm{138b})$

Replacing ${\mu }_{t+}\left(\frac{\nu }{m}\right)$ in Equation 134 with $\xi$, we derive a nonlinear equation in $\xi$ and $p_{\mathrm{\Delta }}^{n+\frac{1}{2}}$:

$G\left(p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right)=-b_1{p}_{\mathrm{\Delta }}^{n+\frac{1}{2}}-c_1\sqrt{\left|p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right|}+b_1{p}_m^{n+\frac{1}{2}}+b_0^{n+\frac{1}{2}}-\frac{S_r}{{\mathrm{\Delta }}_t}\xi =0.(139)$

As in Bilbao et al. (2015b), we have a pair of coupled nonlinear equations, $F(\xi )$, from Equation 67 and $G\left(p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right)$ above. These equations can be solved iteratively for $\xi$ by first solving Equation 67 for $p_{\mathrm{\Delta }}^{n+\frac{1}{2}}$ in terms of $\xi$, substituting the result into Equation 139, and using a Newton-Raphson iteration on the combined equation.

For the linearly implicit quadratized scheme discussed in Section 5.3, we can derive a direct solution for $p_{\mathrm{\Delta }}^{n+\frac{1}{2}}$. Substituting Equation 78 into Equation 139 and using the definition of $\xi$ in Equation 68a, we produce a quadratic equation in $\sqrt{\left|p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right|}$,

$-c_2\text{sign}\left(p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right)\left|p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right|-c_1\text{sign}\left(p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right)\sqrt{\left|p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right|}+c_0=0(140)$

with coefficients dependent on the branching conditions of $g^n$ in Equation 79

$c_2=\left\{\begin{array}{cc}{b}_1+\frac{S_r^2}{a_0},\hfill & h<0\hfill \\ b_1+\frac{S_r^2}{a_2^{n+\frac{1}{2}}{\mathrm{\Delta }}_t},\hfill & h\ge 0\hfill \end{array}\right.(\mathrm{141a})$

$c_1=S_j^{n+\frac{1}{2}}\sqrt{\frac{2}{\rho }},\forall h,(\mathrm{141b})$

$c_0=\left\{\begin{array}{cc}{b}_1{p}_m^{n+\frac{1}{2}}+b_0^{n+\frac{1}{2}}-\frac{S_r}{{\mathrm{\Delta }}_t}\frac{a_1^{n+\frac{1}{2}}}{a_0},\hfill & h<0\hfill \\ b_1{p}_m^{n+\frac{1}{2}}+b_0^{n+\frac{1}{2}}-\frac{S_r}{{\mathrm{\Delta }}_t}\frac{a_1^{n+\frac{1}{2}}-g^n{\mu }_{t+}{\sigma }^n}{a_2^{n+\frac{1}{2}}},\hfill & h\ge 0\hfill \end{array}\right.,(\mathrm{141c})$

with

$a_2^{n+\frac{1}{2}}=a_0+g^{n+\frac{1}{2}}\left(g^{n+\frac{1}{2}}+\frac{{\gamma }_c}{{\mathrm{\Delta }}_t}{\sigma }^n\right).(\mathrm{141d})$

$a_0$ and $a_1^{n+\frac{1}{2}}$ are defined in Equations 68b, c. $c_2,c_1\ge 0$ and we can use the same reasoning as in Harrison-Harsley (2018), Chatziioannou et al. (2019) to guarantee a positive and real solution by

$\text{sign}\left(p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right)=\text{sign}\left(c_0\right),(142)$

and taking the positive solution to the quadratic equation:

$\sqrt{\left|p_{\mathrm{\Delta }}^{n+\frac{1}{2}}\right|}=\frac{-c_1+\sqrt{c_1^2+4c_2\left|c_0\right|}}{2c_2}.(143)$

The value of $p_{\mathrm{\Delta }}^{n+\frac{1}{2}}$ can then be used to compute $\xi$, ${\xi }^{\ast }$, and the associated reed state variables at the next timestep. In the case there is no collision, $g^n$ can be computed as in the top branch of Equation 79 to ensure the contact energy becomes zero. The value of $\xi$ can be used to compute the input flow into the bore using Equation 134.

8.3 Simplified bore synthesis simulation

For our simulation experiment, we designed a simplified bore profile with two toneholes. Through an optimization procedure similar to Lefebvre (2010), the bore length and tonehole position were tuned such that the first impedance peak corresponds to the note C4 when all the holes are closed, and D4 and E4 when one tonehole is open and the other is closed. The resulting bore profile and tonehole positions are shown in Figure 10, and the dots indicate the tonehole positions. The bore radius is defined by the piecewise function

$r\left(z\right)=\left\{\begin{array}{cc}\frac{r_1}{\sqrt{3}},\hfill & 0\ge z<z_1\hfill \\ \tan \left(\theta \right)\left(z_2-z_1\right)+r_1,\hfill & z_1\ge z<z_2\hfill \\ \frac{\beta }{z_a-{\left(z-z_2\right)}^v},\hfill & z_2\ge z<L\hfill \end{array}\right.,(144)$

and is composed of a cylinder, truncated cone, and Bessel horn (Chaigne and Kergomard, 2016). The parameters $\beta$, $z_a$, and $v$ are related to the Bessel horn definition. The radius of the cylinder is chosen to match the missing volume of the truncated cone. The parameters defining the bore as a result of the optimization along with the reed parameters used in the synthesis simulation are given in Table 1.

Figure 10

Figure 10. Simplified bore profile used in final simulation consisting of a cylinder, truncated cone, and a Bessel horn with two toneholes. Tonehole positions are indicated by white markers. Parameter values are provided in Table 1.

Table 1

Table 1. Synthesis parameters used in the final simulation.

The discrete Hamiltonian of the system is the sum of the discrete Hamiltonian from each sub-model

$\mathfrak{h}=+{\mathfrak{h}}_{r,\sigma }+{\displaystyle \sum _{\text{bores}}}{\mathfrak{h}}_b+{\displaystyle \sum _{\text{toneholes}}}{\mathfrak{h}}_{th}.(145)$

Likewise, the dissipated energy is,

$\mathfrak{q}={\mathfrak{q}}_{\sigma }+{\displaystyle \sum _{\text{bores}}}{\mathfrak{q}}_b+{\displaystyle \sum _{\text{toneholes}}}{\mathfrak{q}}_{th}.(146)$

The only boundary term is

$\mathfrak{b}=p_m^{n+\frac{1}{2}}{u}_m^{n+\frac{1}{2}},(147)$

as the boundary terms for each sub-model cancel with each other due to the use of energy-conserving interconnections.

8.4 Results and discussion

The simplified geometry was used to synthesize the first two measures of “Mary Had a Little Lamb,” corresponding to a 4.8 s long audio sample¹. We evaluate the pressure at $l=0$ and times $n+\frac{1}{2}$ corresponding to the boundary condition of the bore. The resulting spectrogram is shown in Figure 11 along with the input mouth pressure, tonehole states, and accumulated energy error. The accumulated energy error is on the order of $10^{-9}$ over the length of the entire simulation. The discrete energy is not an exact sampling of the continuous time energy function due to discretization with the Störmer-Verlet method in the bore model. The gradual increase in accumulated error was determined to be primarily a result of precision loss from computing the synthesis matrices related to the discrete gradient method. Other possible contributions to the shown accumulated energy error drift are the order of operations in the scheme computations and the energy error computation itself (Torin, 2015). The variation of errors is bitwise, displaying only integer multiples of machine precision as in Bilbao and Harrison (2016). It is possible to confirm that the variation in accumulated energy error is random and not systematic via multiple runs of the simulation highlighting the fact that different energy variations are due to the propagation of random finite-precision rounding errors.

Figure 11

Figure 11. Simulation synthesizing the first two measures of “Mary Had a Little Lamb” using the simplified bore described in this section. (A) and (C) display the provided mouth pressure and tonehole states used in the synthesis. (B) and (D) Demonstrate a spectrogram of the resulting synthesized audio and the accumulated energy error.

Simulations were run in Python on an Apple M3 CPU at a 48 kHz sample rate and are not intended for real-time computations. Defining the real-time factor (RTF) as the ratio of the simulated elapsed time and the real processing elapsed time, the simulation performed in this section had an RTF approximately equal to 1/16. An RTF much greater than unity is necessary for real-time computations. A significant bottleneck in the simulation is computation of the viscothermal losses. Running the same simulation without viscothermal losses resulted in an RTF of approximately 1/5.

9 Conclusion

In this article, we have presented a port-Hamiltonian system model of a single-reed woodwind instrument comprised of a lumped reed model, a one-dimensional horn model, and a lumped tonehole model. We then discretized each model using FDTD methods. The PHS framework has aided us in systematically deriving sub-models that can be composed in a modular and energy-conserving fashion, allowing us to construct bores with arbitrary geometries and tonehole placements.

A final simulation based on a simplified bore geometry demonstrates bitwise machine precision variations in the accumulated energy error resulting in a slow long term drift. The drift is not systematic, but random and is based on the propagation of random finite-precision rounding errors. However, further statistical analysis of the scheme should be conducted from hundreds of runs to characterize the behavior of the error as in Hairer et al. (2008). Real-time computation was not a concern in our implementation and proper benchmarking in a compiled programming language, such as C++, would allow for the evaluation of the computational load of the scheme.

Refinements to the lumped reed model and tonehole model were presented. A linearly implicit scheme based on energy quadratization of the Hunt-Crossley contact force was proposed and the behavior of the scheme was compared to a nonlinear implicit scheme solved via a Newton-Raphson iteration. We demonstrated that the linearly implicit scheme maintains the overall characteristic of the interaction. A new low-frequency lumped model of the tonehole was proposed that provides better agreement to existing literature by including damping in the closed tonehole model. The stored energy in the discrete tonehole model is now bounded through use of a switching PHS structure.

The lossy horn equation was discretized using a general Störmer-Verlet scheme, motivating the scheme proposed by Bilbao and Harrison (2016). This initial investigation merits further research into the efficacy of symplectic schemes for discretizing musical systems where numerical dispersion is of a primary concern. A major computational bottleneck in the overall model is due to the complexity of the viscothermal losses as two 16th-order filters are necessary at every spatial index to accurately model the losses (Bilbao et al., 2015a). Development of a more efficient representation remains for future work. Due again to the modular nature of the PHS framework, any future improvements to the sub-models described here or new models of radiation conditions or losses—provided these sub-models are described as a PHS—can be incorporated into the existing model with ease.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

CD: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Writing–original draft, Writing–review and editing. GS: Funding acquisition, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This research was supported by the Vanier Canada Graduate Scholarships program.

Acknowledgments

The authors would like to thank Vasileios Chatziioannou for helpful discussions on symplectic schemes and Mark Rau for discussions regarding the horn equation.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsip.2025.1519450/full#supplementary-material.

Footnotes

¹An audio sample is provided as part of the supplementary material.

References

Avanzini, F., and van Walstijn, M. (2004). Modeling the mechanical response of the reed-mouthpiece-lip system of a clarinet. part i. a one-dimensional distributed model. Acta Acustica united Acustica 90, 537–547.

Google Scholar

Backus, J. (1963). Small-vibration theory of the clarinet. J. Acoust. Soc. Am. 35, 305–313. doi:10.1121/1.1918458

CrossRef Full Text | Google Scholar

Berners, D. P. (1999). Acoustic and signal procesing Techniques for physical Modeling of brass instruments. Ph.D. Thesis. Stanford University.

Google Scholar

Bilbao, S. (2008). “Direct simulation for reed wind instruments,” in 11th International Conference on digital audio effects (DAFx-08). Espoo, Finland.

Google Scholar

Bilbao, S. (2009). Numerical sound synthesis. John Wiley and Sons, Ltd.

CrossRef Full Text | Google Scholar

Bilbao, S., and Chick, J. (2013). Finite difference time domain simulation for the brass instrument bore. J. Acoust. Soc. Am. 134, 3860–3871. doi:10.1121/1.4822479

PubMed Abstract | CrossRef Full Text | Google Scholar

Bilbao, S., and Harrison, R. (2016). Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section. J. Acoust. Soc. Am. 140, 728–740. doi:10.1121/1.4959025

PubMed Abstract | CrossRef Full Text | Google Scholar

Bilbao, S., Harrison, R., Kergomard, J., Lombard, B., and Vergez, C. (2015a). Passive models of viscothermal wave propagation in acoustic tubes. J. Acoust. Soc. Am. 138, 555–558. doi:10.1121/1.4926407

PubMed Abstract | CrossRef Full Text | Google Scholar

Bilbao, S., Torin, A., and Chatziioannou, V. (2015b). Numerical modeling of collisions in musical instruments. Acta Acustica United Acustica, 101.

Google Scholar

Caussé, R., Kergomard, J., and Lurton, X. (1984). Input impedance of brass musical instruments—comparison between experiment and numerical models. J. Acoust. Soc. Am. 75, 241–254. doi:10.1121/1.390402

CrossRef Full Text | Google Scholar

Chaigne, A., and Kergomard, J. (2016). Acoustics of musical instruments. New York, NY: Springer-Verlag.

Google Scholar

Chatziioannou, V. (2010). Forward and inverse modelling of single-reed woodwind instruments with application to digital sound synthesis. Belfast, Northern Ireland: Queen’s University Belfast. Ph.D. thesis.

Google Scholar

Chatziioannou, V. (2019). Structure preserving algorithms for simulation of linearly damped acoustic systems. J. Numer. Analysis, Industrial, Appl. Math. (JNAIAM) 13.

Google Scholar

Chatziioannou, V., Schmutzhard, S., Pàmies-Vilà, M., and Hofmann, A. (2019). Investigating clarinet articulation using a physical model and an artificial blowing machine. Acta Acustica united Acustica 105, 682–694. doi:10.3813/aaa.919348

CrossRef Full Text | Google Scholar

Chatziioannou, V., and van Walstijn, M. (2012). Estimation of clarinet reed parameters by inverse modelling. Acta acustica united Acustica 98, 629–639. doi:10.3813/aaa.918543

CrossRef Full Text | Google Scholar

Dalmont, J.-P., Nederveen, C. J., Dubos, V., Ollivier, S., Méserette, V., and te Sligte, E. (2002). Experimental determination of the equivalent circuit of an open side hole: linear and nonlinear behaviour. Acta Acustica united Acustica 88, 567–575.

Google Scholar

Dubos, V., Kergomard, J., Khettabi, A., Dalmont, J.-P., Keefe, D. H., and Nederveen, C. (1999). Theory of sound propagation in a duct with a branched modal decomposition. Acustica united Acta Acoust. 85, 153–169.

Google Scholar

Ducceschi, M., Bilbao, S., Willemsen, S., and Serafin, S. (2021). Linearly-implicit schemes for collisions in musical acoustics based on energy quadratisation. J. Acoust. Soc. Am. 149, 3502–3516. doi:10.1121/10.0005008

PubMed Abstract | CrossRef Full Text | Google Scholar

V. Duindam, A. Macchelli, S. Stramigioli, and H. Bruyninckx (2009). Modeling and control of complex physical systems: the port-Hamiltonian approach (Springer-Verlag).

Google Scholar

Falaize, A., and Hélie, T. (2016a). Passive guaranteed simulation of analog audio circuits: a port-Hamiltonian approach. Appl. Sci. 273. doi:10.3390/app6100273

CrossRef Full Text | Google Scholar

Falaize, A., and Hélie, T. (2017b). Passive simulation of the nonlinear port-Hamiltonian modeling of a rhodes piano. J. Sound Vib. 390, 289–309. doi:10.1016/j.jsv.2016.11.008

CrossRef Full Text | Google Scholar

Fletcher, N. H. (1993). Autonomous vibration of simple pressure-controlled valves in gas flows. J. Acoust. Soc. Am. 93, 2172–2180. doi:10.1121/1.406857

CrossRef Full Text | Google Scholar

Goldstein, H. (1980). Classical mechanics. 2nd edn. Addison-Wesley.

Google Scholar

Gorrec, Y. L., and Matignon, D. (2013). Coupling between hyperbolic and diffusive systems: a port-Hamiltonian formulation. Eur. J. Control 19, 505–512. doi:10.1016/j.ejcon.2013.09.003

CrossRef Full Text | Google Scholar

Hairer, E., Lubich, C., and Wanner, G. (2000). Geometric Numerical Integration: structure-preserving algorithms for ordinary different equations. Springer.

Google Scholar

Hairer, E., Lubich, C., and Wanner, G. (2003). Geometric numerical integration illustrated by the Störmer–Verlet method. Acta Numer. 12, 399–450. doi:10.1017/s0962492902000144

CrossRef Full Text | Google Scholar

Hairer, E., McLachlan, R. I., and Razakrivony, A. (2008). Achieving brouwer’s law with implicit runge-kutta methods. BIT Numer. Math. 48, 231–243. doi:10.1007/s10543-008-0170-3

CrossRef Full Text | Google Scholar

Harrison-Harsley, R. (2018). Physical modelling of brass instruments using finite-difference time-domain methods. Edinburgh, UK: University of Edinburgh. Ph.D. thesis.

Google Scholar

Hélie, T. (2022). “Elementary tools on Port-Hamiltonian Systems with applications to audio/acoustics,” in 2nd spring School on Theory and Applications of port-Hamiltonian systems (Frauenchiemsee, Germany).

Google Scholar

Hélie, T., and Roze, D. (2016). “Corde non linéaire amortie: fomulation Hamiltonienne à ports, réduction d’ordre exacte et simulation à passivité garantie,” in 13ème Congrès Français d’Acoustique. Le Mans, France.

PubMed Abstract | Google Scholar

Kailath, T. (1980). Linear systems. Prentice-Hall.

Google Scholar

Keefe, D. H. (1982). Theory of the single woodwind tone hole. J. Acoust. Soc. Am. 72, 676–687. doi:10.1121/1.388248

CrossRef Full Text | Google Scholar

Keefe, D. H. (1984). Acoustical wave propagation in cylindrical ducts: Transmission line parameter approximations for isothermal and nonisothermal boundary conditions. J. Acoust. Soc. Am. 75, 58–62. doi:10.1121/1.390300

CrossRef Full Text | Google Scholar

Lefebvre, A. (2010). Computational acoustic methods for the design of woodwind instruments. Montréal, Canada: McGill University. Ph.D. thesis.

Google Scholar

Lefebvre, A., and Scavone, G. P. (2012). Characterization of woodwind instrument toneholes with the finite element method. J. Acoust. Soc. Am. 131, 3153–3163. doi:10.1121/1.3685481

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes, N. (2016). Approche passive pour la modeélisation, la simulation et l’étude d’un banc de test robotisé pour les instruments de type cuivre. Ph.D. thesis, Université Paris 6 (UPMC)

Google Scholar

Lopes, N., and Hélie, T. (2016). Energy balanced model of a jet interacting with a brass player’s lip. Acta Acoust. united Acustica 102, 141–154. doi:10.3813/aaa.918931

CrossRef Full Text | Google Scholar

Lopes, N., Hélie, T., and Falaize, A. (2015). “Explicit second-order accurate method for the passive guaranteed simulation of port-Hamiltonian systems,” in International Federation of Automatic control (IFAC).

Google Scholar

Maschke, B. M., and van der Schaft, A. J. (1993). “Port-controlled Hamiltonian systems: modelling origins and systemtheoretic properties,” in Nonlinear control systems design 1992 (Elsevier), 359–365.

CrossRef Full Text | Google Scholar

McIntyre, M., Schumacher, R., and Woodhouse, J. (1983). On the oscillations of musical instruments. J. Acoust. Soc. Am. 74, 1325–1345. doi:10.1121/1.390157

CrossRef Full Text | Google Scholar

Mignot, R., Hélie, T., and Matignon, D. (2008). “Stable relization of a delay system modeling a convergent acoustic cone,” in 16th Mediterranean Conference on Control and Automation. Ajaccio, France.

Google Scholar

Mignot, R., Hélie, T., and Matignon, D. (2010). Digital waveguide modeling for wind instruments: Building a state–space representation based on the Webster–Lokshin model. IEEE Trans. Audio, Speech, Lang. Process. 18, 843–854. doi:10.1109/tasl.2009.2038671

CrossRef Full Text | Google Scholar

Müller, R. (2021). Time-continuous power-balanced simulation of nonlinear audio circuits: realtime processing framework and aliasing rejection. Sorbonne Université: Ph.D. thesis.

Google Scholar

Müller, R., and Hélie, T. (2021). “Fully-implicit algebro-differential parameterization of circuits,” in Proceedings of the 23^rdInternational Conference on digital audio effect (DAFx2020).

Google Scholar

Sanz-Serna, J., and Calvo, M. (1994). Numerical Hamiltonian problems. London, UK: Chapman and Hall.

Google Scholar

Scavone, G. (1997). An acoustic analysis of single-reed woodwind instruments with an Emphasis on design and performance issues and digital waveguide modeling techniques. Stanford, CA: Stanford University. Ph.D. thesis.

Google Scholar

Scavone, G., and Cook, P. (1998). “Real-time computer modeling of woodwind instruments,” in 1998 International Symposium on musical acoustics Leavenworth, WA.

Google Scholar

Scavone, G. P. (2024). “An open-source project for wind instrument modeling using digital waveguides,” in 186^thMeeting of the acoustical Society of America.

Google Scholar

Shen, J., Xu, J., and Yang, J. (2018). The scalar auxiliary variable (SAV) approach for gradient flows. J. Comput. Phys. 353, 407–416. doi:10.1016/j.jcp.2017.10.021

CrossRef Full Text | Google Scholar

Smith, J. O. (1986). “Efficient simulation of the reed-bore and bow-string mechanisms,” in Proceedings of the 12th International Conference on music computing (Netherlands: Den Haag) Available at: http://books.w3k.org/

Google Scholar

Smith, J. O. (1992). Physical modeling using digital waveguides. Comput. Music J. 16, 74–91. doi:10.2307/3680470

CrossRef Full Text | Google Scholar

Smith, J. O. (2010). Physical audio signal processing. W3K Publishing.

Google Scholar

Torin, A. (2015). Percussion instrument modelling in 3D: sound synthesis through time domain numerical simulation. Edinburgh, UK: University of Edinburgh. Ph.D. thesis.

Google Scholar

Välimäki, V. (1995). Discrete-time modeling of acoustic tubes using fractional delay filters. Helsinki, Finland: Helsinki University of Technology. Ph.D. thesis.

Google Scholar

van der Schaft, A. (2006). “Port-Hamiltonian systems: an introductory survey,” in Proceedings of the International Congress of Mathematicians. Madrid, Spain.

Google Scholar

van der Schaft, A., and Jeltsema, D. (2014). Port-Hamiltonian systems theory: an introductory Overview. Hamilt. Syst. Theory Introd. Overv. 1, 173–378. doi:10.1561/2600000002

CrossRef Full Text | Google Scholar

van Walstijn, M., and Avanzini, F. (2007). Modelling the mechanical response of the reed-mouthpiece-lip system of a clarinet. part ii: a lumped model approximation. Acta Acustica united Acustica 93, 435–446.

Google Scholar

van Walstijn, M., Bridges, J., and Mehes, S. (2016). “A real-time synthesis oriented tanpura model,” in Proceedings of the 19^thInternational Conference on digital audio effects (DAFx-16).

Google Scholar

van Walstijn, M., and Campbell, M. (2003). Discrete-time modeling of woodwind instrument bores using wave variables. J. Acoust. Soc. Am. 113, 575–585. doi:10.1121/1.1515776

PubMed Abstract | CrossRef Full Text | Google Scholar

van Walstijn, M., Chatziioannou, V., and Athanasopoulos, N. (2024a). An explicit scheme for energy-stable simulation of mass-barrier collisions with contact damping and dry friction. IFAC-PapersOnLine 58, 214–219. doi:10.1016/j.ifacol.2024.08.283

CrossRef Full Text | Google Scholar

van Walstijn, M., Chatziioannou, V., and Bhanuprakash, A. (2024b). Implicit and explicit schemes for energy-stable simulation of string vibrations with collisions: Refinement, analysis and comparison. J. Sound Vib. 569, 117968. doi:10.1016/j.jsv.2023.117968

CrossRef Full Text | Google Scholar

van Walstijn, M., and Scavone, G. (2000). “The wave digital tonehole model,” in 2000 International computer music Conference Berlin, Germany.

Google Scholar

Wendlandt, J. M., and Marsden, J. E. (1997). Mechanical integrators derived from a discrete variational principle. Phys. D. 106, 223–246. doi:10.1016/s0167-2789(97)00051-1

CrossRef Full Text | Google Scholar

Wetzel, V., Hélie, T., and Silve, F. (2019). “Power balanced time-varying lumped parameter model of a vocal tract: modeling and simulation,” in 26th International Conference on sound and Vibration.

Google Scholar

Willemsen, S. (2021). The Emulated ensemble: real-time simulation of musical instruments using finite-difference time-domain methods. Aalborg, Denmark: Aalborg University. Ph.D. thesis.

Google Scholar

Yalçin, Y., Sümer, L. G., and Kurtulan, S. (2015). Discrete-time modeling of Hamiltonian systems. Turkish J. Electr. Eng. Comput. Sci. 23, 149–170. doi:10.3906/elk-1212-23

CrossRef Full Text | Google Scholar

Yang, X. (2016). Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends. J. Comput. Phys. 327, 294–316. doi:10.1016/j.jcp.2016.09.029

CrossRef Full Text | Google Scholar

Zhao, J., Want, Q., and Yang, X. (2016). Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach. Int. J. Numer. Methods Eng. 110, 279–300. doi:10.1002/nme.5372

CrossRef Full Text | Google Scholar

Zwikker, C., and Kosten, C. W. (1949). Sound absorbing materials. New York: Elsevier.

Google Scholar