Physical Modelling of Piano Sound
Haifan Xie

UTF8


Physical Modelling of Piano Sound
    Haifan XieThird-year master student, School of Accounting, Guangdong University
of Foreign Studies, Guangzhou, China. Email: fan455@outlook.com.
    Received XXX; accepted ZZZ
====================================================================================================================================================

§ ABSTRACT 
This paper develops a comprehensive physical model and numerical implementation
schemes for a grand piano, building upon the prior works of Chabassier
et al. The model encompasses various subsystems, including hammer
felt, hammer shank, string, soundboard, air and room barriers, each
modelled in three dimensions to approach their realistic dynamics.
A general framework for 3D elastic solids accounting for prestress
and prestrain is introduced, particularly addressing the complexities
of prestressed piano strings. The study also examines coupling between
subsystem through mechanisms of surface force transmission and displacement/velocity
continuity. To facilitate numerical simulations, strong PDEs are translated
into weak ODEs via a flexible space discretization approach. Modal
transformation of system ODEs is then employed to decouple and reduce
DOFs, and an explicit time discretization scheme is customized for
generating digital audio in the time domain. The study concludes with
a discussion of the piano model’s capabilities, limitations, and potential
future enhancements.




§ INTRODUCTION


This study aims to present a detailed physical model of a grand piano
and its numeric implementation schemes. The main features of our model
and numeric schemes include:


  * As 3D as possible. Subsystems including hammer felt, hammer shank,
soundboard, air and room barriers are all fully 3D geometrically modelled.
This also means unless for rigid bodies (hammer shank), their 3-directional
displacements fully vary with their 3D positions. Our geometric model
of piano strings is semi-3D, viz. a cylinder geometry described by
coordinates of one axis, but accounts for the strings' 3-directional
displacements. 

  * A general framework for 3D elastic solid. Rooted in the 3D elasticity
theory, this framework not only follows the classic stress-strain
relationship, but also introduces a simple yet intuitive model for
prestrain and prestress. For significantly prestressed structures
like piano strings, we offer a more “naturally linear” prestress
model compared with the elaborated geometrically exact nonlinear stiff
string model in <cit.>.

  * Consideration of coupling between subsystems of the piano. Two mechanisms
are considered here: surface force transmission and displacement/velocity
continuity. Furthermore, nonlinear coupling due to collision rather
than fixation is considered in the interaction between hammer felt
and string, and discussed (but failed to implement) in the interaction
between string and two bridge pins.

  * A straightford framework for weak forms and system ordinary differential
equations (ODE). For the sake of numeric simulation via the finite
element method (FEM), we transform the strong form partial differential
equations (PDE) into first or second order system ODEs for each subsystem
of the piano. This is achieved via a flexible space discretization
framework, which can easily incoporate Dirichlet boundary conditions.

  * Modal transformation for solving system ODEs. By solving generalized
eigenvalue problems of the mass and stiffness matrices, we decouple
the coupled system ODEs and significantly reduce the large vector
of degrees of freedom (DOF) into a much smaller vector of modal DOFs.

  * An explicit time discretization scheme customized for solving modal
ODEs exhibiting source term coupling. This scheme is relatively efficient
in that no mid steps are evaluated and no inversion of large non-diagonal
matrices is needed. This scheme is relatively accurate in that rhs
source terms of the next time step are used for approximate integration
whenever possible, and more than one iteration of numeric integration
at each time step may be run to improve convergence.

The rest of this paper starts from a basic 3D elastic prestressed
solid model in section 2, which will be applied in several cases later.
Subsequently sections 3 to 5 introduce models for subsystems of the
piano: soundboard, string, air and room barriers. Section 6 investigates
the mechanisms for coupling between these subsystems, where hammer
felt and shank models are introduced. Section 7 presents the modal
superposition method and explicit time stepping schemes to solve the
derived system ODEs numerically, and a comprehensive computation procedure
for running the simulation. Section 8 gives a brief summary of the
whole model and discusses current limitations and future outlooks.

Here we explain some notation patterns that will appear throughout
this paper. If not explicitly specified, symbols defined only apply
to the current section. In cases where avoiding symbol conflict is
necessary, superscripts ^(a),^(b),^(c),^(d),^(e),^(f)
refer to the string, soundboard, air, hammer shank, hammer felt parts,
room barriers of a piano physical system respectively. For the convenience
of notation, (x,y,z) and (x_1,x_2,x_3), (u,v,w) and
(u_1,u_2,u_3) are used interchangeably.



§ PRELIMINARY: 3D LINEAR ELASTIC SOLID WITH PRESTRESS




In this study, the physical models of piano strings and soundboard
are based on the full or reduced versions of the linear theory of
elasticity <cit.>. In a 3D Cartesian coordinate
system (x,y,z), consider a material defined over a space Ω⊂ℝ^3
with a boundary Γ⊂ℝ^3, with homogenous or
heterogenous density ρ(x,y,z), and dynamic displacements in
the 3 directions as 

    u(x,y,z,t)=[u(x,y,z,t),v(x,y,z,t),w(x,y,z,t)]^⊤.

Based on the displacement field, we shall perform a force analysis
of the material, considering two kinds of forces: surface force and
body force.

Stress is a major source of surface force for elastic material, and
can be derived from the stress-strain relation. The strain matrix
(second-order symmetric tensor) is expressed by the gradients of displacements
as 

    ϵ=[[ ϵ_11 ϵ_12 ϵ_13; ϵ_12 ϵ_22 ϵ_23; ϵ_13 ϵ_23 ϵ_33 ]]=[[           ∂_xu 1/2(∂_yu+∂_xv) 1/2(∂_zu+∂_xw); 1/2(∂_yu+∂_xv)           ∂_yv 1/2(∂_zv+∂_yw); 1/2(∂_zu+∂_xw) 1/2(∂_zv+∂_yw)           ∂_zw ]],

which is the linearized Green-Lagrange strain. It can also be written
in vector form as 

    ϵ=[[  ϵ_11;  ϵ_22;  ϵ_33; 2ϵ_12; 2ϵ_13; 2ϵ_23 ]]=[[      ∂_xu;      ∂_yv;      ∂_zw; ∂_yu+∂_xv; ∂_zu+∂_xw; ∂_zv+∂_yw ]]=∑_i=1^3H_i∇ u_i,

where 

    H_1=[[ 1 0 0; 0 0 0; 0 0 0; 0 1 0; 0 0 1; 0 0 0 ]], H_2=[[ 0 0 0; 0 1 0; 0 0 0; 1 0 0; 0 0 0; 0 0 1 ]], H_3=[[ 0 0 0; 0 0 0; 0 0 1; 0 0 0; 1 0 0; 0 1 0 ]].

From the generalized Hooke's law, the stress-strain relationship is
(here the stress is Cauchy stress) 

    σ=[[ σ_11; σ_22; σ_33; σ_12; σ_13; σ_23 ]]=Dϵ, D=[[ D_11 D_12 D_13 D_14 D_15 D_16;      D_22 D_23 D_24 D_25 D_26;           D_33 D_34 D_35 D_36;                D_44 D_45 D_46;                     D_55 D_56; Sym.                     D_66 ]],

and the stress matrix (second-order symmetric tensor) is defined as

    σ=[[ σ_11 σ_12 σ_13; σ_12 σ_22 σ_23; σ_13 σ_23 σ_33 ]].

As shown in figure <ref>, there are 9 stress forces
acting on the 3 positive surfaces (as well as 9 stress forces on the
3 negative surfaces not shown) of an infinitesimal volume of material.
D is called the constitutive matrix and its inverse
D^-1 is called the compliance matrix. The i th
row or column of the stress matrix can be expressed as 

    σ_i=∑_j=1^3A_ij∇ u_j, A_ij=H_i^⊤DH_j, A=[[ A_1; A_2; A_3 ]]=[[ A_11 A_12 A_13; A_21 A_22 A_23; A_31 A_32 A_33 ]],

where block matrix A is symmetric. As σ_i
represents surface forces, converting it to body force would result
in ∇·σ_i. 

Prestress is considered as the initial stress of material in static
equilibrium, which is not included in the above analysis of stress.
Incorporating prestress into our analysis would require some additional
work as presented in appendix <ref>,
which we shall refer to. Define the static tension field matrix (second-order
symmetric tensor) and its vector form as 

    T(x,y,z)=[[ T_11 T_12 T_13; T_12 T_22 T_23; T_13 T_23 T_33 ]], T(x,y,z)=[[ T_11; T_22; T_23; T_12; T_13; T_23 ]],

which can be visualized by figure <ref> similar
to stress. Since T is the prestress in static equilibrium,
∇·T=0 should hold (ignoring
gravity). As per (<ref>), the contribution of
prestress to the total stress consists of the static prestress T,
and the dynamic prestress τ which is also a symmetric
tensor. Denote T_i and τ_i
as the i th row or column of T and τ
respectively, then we find 

    τ_i=∑_j=1^3B_ij∇ u_j, B_ij=H_i^⊤DΨH_j, B=[[ B_1; B_2; B_3 ]]=[[ B_11 B_12 B_13; B_21 B_22 B_23; B_31 B_32 B_33 ]],

where block matrix B is non-symmetric. Similar to
σ_i previously discussed, T_i
and τ_i contain the 3 sources of prestress in
the same i th axis. For static equilibrium, ∇·T_i=0
should hold; then in dynamic states, ∇·τ_i
is the body force of prestress in the i th axis.

Besides stress and prestress which are conservative, non-conservative
forces like damping force often exist. For elastic non-metallic solid,
viscoelastic damping is often the predominant damping <cit.>.
We hereby adopt a simple viscoelastic model: the viscous damping force
is positively proportional to the first-order time derivative of stress
and prestress[An explanation for this is that stress and prestress, along with viscous
damping force, are resistance to deformation, and so their relations
should be positive. ]. Occurring on surfaces of small volumes, the damping forces can be
visualized by figure <ref> similar to stress. The
vector of damping forces parallel to the x_i axis is defined
as 

    ς_i=2μ∂_t(σ_i+τ_i),

where 2μ>0 is the damping coefficient. Other damping models can
also be incoporated into our model, e.g. structural damping that adds
frequency-dependent imaginary parts to elasticity coefficients <cit.>.
As for non-conservative forces besides damping, we represent them
as a single force F=[F_1,F_2,F_3]^⊤ and
will investigate them for different physical systems later.

Having derived all the relevant forces acting on an infinitesimal
volume of solid material, we can invoke Newton's second law to derive
3 PDEs as 

    ρ∂_ttu_i=∇·G_i+F_i, i=1,2,3.

where vector G_i=σ_i+τ_i+ς_i
is defined as the i th row or column of symmetric tensor G
representing all surface forces except the static prestress. The above
equation can also be derived via a Lagrange formulation of the virtual
work and the variations of kinetic and potential energy. We then seek
for a weak form of it via variational formulation. Multiplying an
arbitrary test function vector ψ(x,y,z) on both
sides of (<ref>) and integrating over Ω yields 

    ∫_Ωψρ∂_ttu_idV-∫_Ωψ∇·G_idV     =∫_ΩψF_idV
    ∫_Ωψρ∂_ttu_idV+∫_Ω(∇ψ)G_idV     =∫_ΩψF_idV

where dV=dxdydz and applying
the gradient operator ∇ to a vector results in its Jacobian
matrix. In the above formulation, integration by parts is utilized
with Neumann boundary condition imposed as 

    G_i·dΓ=0,

where dΓ is the outward normal vector
of the tangent plane of any point on Γ. Note that this condition
applies only when the test function is non-zero on the boundary, viz.
Neumann boundary conditions need not be satisfied at points where
Dirichlet boundary conditions are present. 

In the next following sections, we shall apply the above 3D elastic
solid material model to various parts of the piano physical system,
including soundboard, string, room material and hammer felt. Also,
we will use space discretization and the Galerkin method to transform
(<ref>) into weak forms of second-order ODE.



§ MODEL FOR PIANO SOUNDBOARD




The symbols in this section follow section 2 if not explicitly defined.
Definitions for scalars, vectors and matrices with subscript _0,
if not explicited stated, are automatically inferred from their counterparts
without subscript _0. 

The study treats the grand piano soundboard as a multi-layer plate
with irregular geometry, and physically models it as a 3D structure.
In <cit.>,
the soundboard was modelled as a Reissner-Mindlin plate, accounting
for the ribs and bridges by making the thickness, density, and elastic
coefficients position-dependent; however, it is unclear how the different
orthotropic angles of each layer were handled. In <cit.>,
the soundboard was modelled as a Kirchhoff-Love plate, considering
a 90 degrees orthotropic rotation of the ribs; however, as the Kirchhoff-Love
plate specifies only 1 unknown for the 3D displacement field, it may
not be accurate enough, especially regarding the ribs and bridges
which are much thicker than the board. Moreover, the role of rim and
lid in shaping the piano timbre seems under-evaluated.

To better account for the soundboard's multi-layer feature, the current
study implements a fully 3D soundboard model. Figure <ref>
shows a simplified soundboard sketch based on <cit.>.
In our physical model, the soundboard is an entity composed of 5 parts:
board, ribs, bridges, rim, lid. The board is parallel to the xOy
plane, with the ribs attached to its under side and the bridges attached
to its upper side. The horizontal fibers of the board and bridges
are approximately parallel, while the horizontal fibers of ribs are
almost orthogonal to that of the board. The rim consists of 2 parts:
the inner rim connects the board's boundary from below using bolts
and dowels, and the outer rim encases the board and the lower rim.
The lid, as large as the board and with an angle to it, is an extension
of the outer rim through a wooden stick (considered part of the lid)
and hinges. All parts of the soundboard are treated as a whole in
computation, meaning that the displacements of part intersections
are uniform.

We define the occupied space of soundboard as Ω⊂ℝ^3,
and the boundaries of soundboard as Γ_1,Γ_2⊂Γ⊂ℝ^3.
Here Γ contains all the 2D surfaces of the 3D volume of soundboard;
Γ_1 contains only the underside surfaces of the inner and
outer rims, as marked red at the right of figure 2; Γ_2
is Γ excluding Γ_1, the vibrating parts. We shall
see in the following that Γ_1 and Γ_2 are where
the Dirichlet boundary conditions apply to the soundboard and the
air respectively.

Based on (<ref>), define the displacement field of
the soundboard as u^*≈u=[u,v,w]^⊤
(here superscript ^* means the exact solution) as 

    u(x,y,z,t)     =φ_1(x,y,z)·[S_1ξ(t)+S_0,1ξ_0(t)],
    
    v(x,y,z,t)     =φ_2(x,y,z)·[S_2ξ(t)+S_0,2ξ_0(t)],
    
    w(x,y,z,t)     =φ_3(x,y,z)·[S_3ξ(t)+S_0,3ξ_0(t)],

where the sizes of vectors and matrices are 

    ξ:N×1, ξ_0:N_0×1,
       φ_1:N_1×1, S_1:N_1× N, S_0,1:N_0,1× N_0,
       φ_2:N_2×1, S_2:N_2× N, S_0,2:N_0,2× N_0,
       φ_3:N_3×1, S_3:N_3× N, S_0,3:N_0,3× N_0.

In (<ref>), each scalar displacement unknown
is expressed as the dot product of a space function vector and time
function vector. In FEM, φ_i is the vector
of shape functions (also known as interpolation functions) for variable
u_i, and ϑ_i=S_iξ+S_0,iξ_0
is the vector of coefficients of shape functions, often interpreted
as nodal displacements. The only unknown here is ξ(t),
the vector of unknown DOFs, whereas other vectors and matrices are
known. ξ(t) is defined so that through some linear
transformation by S_1 and S_0,1ξ_0(t),
the nodal displacements ϑ_i can be obtained.
The reason we do not simply define u_i=φ_i·ξ_i
but consider a linear transformation is that it provides additional
flexibility when imposing Dirichlet boundary conditions. For instance,
when u+v rather than u or v is known to be some nonzero functions
on the space boundary, ξ only needs to contain u
and u+v is incorporated in ξ_0. We can thus
see that total DOF, viz. the number of time functions we need to solve,
is N that does not necessarily equal to N_1+N_2+N_3. Often
in simple cases, matrix S=[S_1^⊤,S_2^⊤,S_3^⊤]^⊤
is diagonal with many ones and some zeros on the diagonal.

To express the displacement and its gradient in the DOF vector, we
find the below relations:

    u_i=P_iξ+P_0,iξ_0, ∇ u_i=Q_iξ+Q_0,iξ_0,
       P=[[ P_1; P_2; P_3 ]]=[[ φ_1^⊤S_1; φ_2^⊤S_2; φ_3^⊤S_3 ]], P_0=[[ P_0,1; P_0,2; P_0,3 ]]=[[ φ_1^⊤S_0,1; φ_2^⊤S_0,2; φ_3^⊤S_0,3 ]],
       Q=[[ Q_1; Q_2; Q_3 ]]=[[ (∇φ_1)^⊤S_1; (∇φ_2)^⊤S_2; (∇φ_3)^⊤S_3 ]], Q_0=[[ Q_0,1; Q_0,2; Q_0,3 ]]=[[ (∇φ_1)^⊤S_0,1; (∇φ_2)^⊤S_0,2; (∇φ_3)^⊤S_0,3 ]].

Given that all parts of the soundboard are mainly made of wood, an
orthotropic material whose orthotropic directions are determined by
fibers <cit.>, the compliance matrix writes 

    D_orth^-1=[[     1/E_x -ν_xy/E_x -ν_xz/E_x         0         0         0;               1/E_y -ν_yz/E_y         0         0         0;                         1/E_z         0         0         0;                                    G_xy         0         0;                                              G_xz         0;      Sym.                                              G_yz ]].

Note that different layers would have different elastic coefficients.
It is then necessary to consider that for a certain layer like ribs,
the global axes x,y may need to be rotated by an angle α,
as denoted in figure <ref>, to become the material
orthotropic axes x',y' and for a certain layer. The actual constitutive
matrix is then 

    D=Z^⊤D_orthZ, Z=[[     C^2     S^2       0      SC       0       0;     S^2     C^2       0     -SC       0       0;       0       0       1       0       0       0;    -2SC     2SC       0 C^2-S^2       0       0;       0       0       0       0       C       S;       0       0       0       0      -S       C ]],

where S=sinα,C=cosα. Readers may refer to <cit.>
for a detailed deduction. After rotation, the constitutive matrix
would have 13 non-zero entries for its upper triangular part. 

Combining (<ref>) with (<ref>)
(<ref>) , the following relations are found:

    σ_i    =A_iQξ+A_iQ_0ξ_0,
    τ_i    =B_iQξ+B_iQ_0ξ_0,

Now we subsitute (<ref>) (<ref>)
into (<ref>) (<ref>) (<ref>),
and use P_i^⊤ as a vector of test functions
for the i th PDE. These would yield 3 groups of weak-form equations,
each group having N equations. Summing the 3 groups into 1 group,
we derive a system of N coupled second-order ODEs as 

    Mξ̈(t)+Cξ̇(t)+Kξ(t)=f(t),

where 

    M=∫_ΩρP^⊤PdV, M_0=∫_ΩρP^⊤P_0dV,
       K=∫_ΩQ^⊤(A+B)QdV, K_0=∫_ΩQ^⊤(A+B)Q_0dV,
       C=2μK, C_0=2μK_0,
       f(t)=∫_ΩP^⊤FdV-M_0ξ̈_0(t)-C_0ξ̇_̇0̇(t)-K_0ξ_0(t),

are called the mass matrix, stiffness matrix, damping matrix and force
vector respectively. In order that reasonable solutions can be found,
equation (<ref>) is constrained by the following Dirichlet
boundary condition and initial conditions: 

    u|_x∈Γ_1=0,


    u|_t=0=∂_tu|_t=0=0.

The choice of this boundary condition stems from that the whole soundboard
is hard supported by beams and legs on the under side of rim (see
figure 2). Thanks to the flexiblility of our 3D model, the rim can
be treated as a natural extension of the main board, thereby capturing
more nuanced boundary conditions compared to the clamped or simply
supported cases in 2D models. Given approriate boundary conditions,
the mass and stiffness matrices are symmetric positive definite.



§ MODEL FOR PIANO STRINGS




The symbols in this section follow sections 2 and 3 if not explicitly
defined. Definitions for scalars, vectors and matrices with subscript
_0, if not explicited stated, are automatically inferred from
their counterparts without subscript _0. 

The study treats a grand piano string as a cylinder material, and
physically models it as a 1D beam prestressed longitudinally. As shown
in figure <ref>, the string is fixed at the agraffe end
(immobile) and coupled to the bridge end (mobile) <cit.>
through pins. Define that the central line of string range start from
the tunning pin [-L_0,0,0]^⊤, then go through the agraffe
x_0=[0,0,0]^⊤, the front bridge pin x_1=[L_1,0,0]^⊤,
the rear bridge pin x_2=[L_2,0,0]^⊤, and
finally the hitch pin x_3=[L_3,0,0]^⊤, where
L_1 is the so-called “speaking length”; let r be the radius
of cross-section and ρ be the homogenous density per unit volume.
It is known that the vibration of piano strings is primarily vertical
(z direction) <cit.>. Nevertheless, logitudinal
vibration (x direction) <cit.> and horizontal
vibration (y direction) <cit.> may also be essential,
as they contribute respectively to the sound precursor <cit.>
and double-decay <cit.> phenomena. Two more kinds of
vibration that require attention are the rotations of the string's
cross-section towards the y and z axes, so as to account for
the stiffness of piano strings that may result in slightly inharmonic
sounds <cit.>. To incorporate all these essential
kinds of vibrations, define the displacement field of a piano string
as u^*≈u=[u,v,w]^⊤ in a
similar notion to (<ref>) (thus some details
are omitted here) as 

    u(x,y,z,t)     =φ_1(x)·[S_1ξ(t)+S_0,1ξ_0(t)]-yφ_4(x)·[S_4ξ(t)+S_0,4ξ_0(t)]
        -zφ_5(x)·[S_5ξ(t)+S_0,5ξ_0(t)],
    
    v(x,y,z,t)     =φ_2(x)·[S_2ξ(t)+S_0,2ξ_0(t)],
    
    w(x,y,z,t)     =φ_3(x)·[S_3ξ(t)+S_0,3ξ_0(t)],

where subscripts _1,_2,_3,_4,_5 refer to vectors or matrices
defined for longitudinal, horizontal, vertical, horizontal rotational
and vertical rotational vibrations respectively, as marked u,v,w,α,β
in figure <ref>. From a 3D perspective, it can be seen
that with regard to the y and z axes, the x displacement
is first-order modelled, while the y and z displacements are
zeroth-order modelled. For the rest of this section, we only discuss
the particularities of piano string model relative to those presented
in sections 2 and 3. 

The below relations hold 

    u_i    =P_iξ+P_0,iξ_0, ∇ u_i=Q_iξ+Q_0,iξ_0,
    P    =[[ P_1; P_2; P_3 ]]=[[ φ_1^⊤S_1-yφ_4^⊤S_4-zφ_5^⊤S_5;                     φ_2^⊤S_2;                     φ_3^⊤S_3 ]],
    Q    =[[ Q_1; Q_2; Q_3 ]]=[[ (∇φ_1)^⊤S_1-y(∇φ_4)^⊤S_4-z(∇φ_5)^⊤S_5;                           (∇φ_2)^⊤S_2;                           (∇φ_3)^⊤S_3 ]],

where P_0 and Q_0 are defined
similar to P and Q. From (<ref>)
(<ref>), it can be found that ϵ_22=ϵ_33=2ϵ_23=0,
thus the strain vector reduces to ϵ=[ϵ_11,2ϵ_12,2ϵ_13]^⊤,
and (<ref>) can be redefined as 

    H_1=[[ 1 0 0; 0 1 0; 0 0 1 ]], H_2=[[ 0 0 0; 1 0 0; 0 0 0 ]], H_3=[[ 0 0 0; 0 0 0; 1 0 0 ]].

As the steel used to manufacture piano strings can be considered as
isotropic material, the stress-strain relation in (<ref>)
reduces to 

    σ=[[ σ_11; σ_12; σ_13 ]]=[[ D 0 0; 0 G 0; 0 0 G ]][[  ϵ_11; 2ϵ_12; 2ϵ_13 ]]=Dϵ,

where D=E(1-ν)/(1+ν)(1-2ν), G=E/2(1+ν);
E is the young's modulus, ν is the Poisson's ratio, and G
is the shear modulus. 

It now suffices to derive a system of N coupled second-order ODEs
simply by substituting (<ref>) into (<ref>)
(<ref>). In order that reasonable solutions
can be found, equation (<ref>) is at least constrained by
the following Dirichlet boundary conditions: 

    u|_x=0=u|_x=L_3=0,

which means the motion of string vanishes at the agraffe point and
the hitch point. Another Dirichlet boundary condition at the coupling
point x_1 will be discussed in section <ref>.
The initial conditions are

    u|_t=0=∂_tu|_t=0=0.





§ MODEL FOR SOUND RADIATION IN THE AIR




 §.§ Model for the air




The study adopts a 3D acoustic wave model to simulate piano sound
radiation in the air. For simulating wave propagation from the piano
soundboard to a listener situated at a particular location, <cit.>
used Rayleigh integral to compute the acoustic pressure field, which
is able to simulate the different delays and decays of sound at different
positions, but unable to account for boundary conditions. <cit.>
used 2 first-order linear acoustic equations for 4 coupled unknowns:
acoustic pressure (1 unknown) and acoustic velocity (3 unknowns).
This approach excels at capturing reflections at spatial boundaries
like walls and soundboard-air coupling, but seems unable to produce
decaying room impulse response due to the absence of damping terms[It seems the acoustic model in <cit.> relies on
soundboard-air coupling to produce the damping. When the acoustic
equations are fully coupled to the soundboard equations, the soundboard's
damping mechanisms may apply to the air as well. This means even if
the acoustic equations contains no damping, the overall energy may
be stll stable and waves with infinite amplitudes would not occur.
However, we choose to add viscous damping terms to the acoustic equations
for two reasons. Firstly, it can be observed in audio recording that
a very short impulse in a room gets a decaying response, which attributes
to the reverberation effect. Secondly, the coupling computation approach
we shall introduce later actually relies on decoupling and iterative
strategies, meaning that without the room itself's damping mechanisms
the computed pressure may exhibit infinite amplitudes.]. 

Fluid dynamics exhibit both linear and nonlinear behaviours as decribed
by the Navier-Stokes equation, but considering that sound pressure
fluctuations in air are often small, linearized models for sound radiation
in the air can achieve a satisfactory level of accuracy for us. The
acoustic radiation equations we utilize are based on the linearized
fluid dynamics equations <cit.>, including conservation
of mass and momentum. It works for isotropic, compressible, viscous,
adiabatic flow with homogenous ambient states, in Eulerian description.
Damping due to heat conduction is simplified away here and readers
may refer to <cit.> for incorporating
thermal damping. The governing PDEs write 

    1/ρ c^2ṗ+∇·u=0,
       ρu̇=∇·G+F,

where u(x,y,z,t) is the acoustic velocity field; p(x,y,z,t)
is the perturbation of acoustic pressure field, also should be the
final digital audio signal; ρ is the air density; c is the
sound propagation speed in the air; G is the symmetric
surface force tensor; the gradient operator ∇ applied to a
vector returns its Jacobian matrix, and the divergence operator ∇·
applied to a matrix returns a vector of each row's divergence; F=[F_1,F_2,F_3]^⊤
is the external force (viz. body force) treated as zero in this study.
According to the Navier-Stokes equation, the surface force tensor
G is defined as 

    G=-pI+μ_1(∇u̇+(∇u̇)^⊤)+μ_2(∇·u̇)I,

where μ_1 and μ_B are dynamic viscosity and bulk viscosity
coefficients accounting for energy dissipation, and μ_2=μ_B-2/3μ_1.
To express the divergence of surface force using acoustic velocity
gradients, we find 

    ∇·G_i=-∇ p+∇·ς_i, ς_i=∑_j=1^3(μ_1A_ij+μ_2B_ij)∇ u_j,
       A=[[ A_1; A_2; A_3 ]]=[[ A_11 A_12 A_13; A_21 A_22 A_23; A_31 A_32 A_33 ]], B=[[ B_1; B_2; B_3 ]]=[[ B_11 B_12 B_13; B_21 B_22 B_23; B_31 B_32 B_33 ]],
       A=[[ 2 0 0 0 0 0 0 0 0; 0 1 0 1 0 0 0 0 0; 0 0 1 0 0 0 1 0 0; 0 1 0 1 0 0 0 0 0; 0 0 0 0 2 0 0 0 0; 0 0 0 0 0 1 0 1 0; 0 0 1 0 0 0 1 0 0; 0 0 0 0 0 1 0 1 0; 0 0 0 0 0 0 0 0 2 ]], B=[[ 1 0 0 0 1 0 0 0 1; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 1 0 0 0 1 0 0 0 1; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 1 0 0 0 1 0 0 0 1 ]],

where G_i, ς_i is the
i th column or row of G, ς;
here ς is the damping force tensor. The 2
PDEs (<ref>) (<ref>) then become[It is also viable to decouple the acoustic equations into a single
second-order equation containing only u as unknown.
We choose to not do so because it would result in a damping matrix
not diagonalizable by the mass or stiffness matrices. To fully decouple
the system ODEs would then require doubling the DOFs (equivalent to
6 variables) to transform into a first-order system. This appears
suboptimal to us as the original first-order system of acoustic equations
has only 4 varaibles.] 

    ρu̇_i-∇·ς_i+∂_x_ip=F_i (i=1,2,3),
       ṗ+ρ c^2∇·u=0,

where the divergence of velocity can be expressed in velocity gradients
as 

    ∇·u=∑_j=1^3H_j∇ u_j,
       H=[[ H_1 H_2 H_3 ]]=[[ 1 0 0 0 0 0 0 0 0; 0 0 0 0 1 0 0 0 0; 0 0 0 0 0 0 0 0 1 ]].

For acoustic modes the vorticity ∇×u=0
may be assumed zero so that a decoupled equation containing only p
as unknown may be obtained, but we choose to not do so because the
velocity field of any solid to couple with may not satisfy zero vorticity.

Now we seek for the weak forms of (<ref>)
(<ref>) via the Galerkin method. Multiplying
an arbitrary test function vector ψ on both sides
of (<ref>) yields 

    ∫_Ωψρu̇_idV-∫_Ωψ(∇·ς_i)dV     +∫_Ωψ∂_x_ipdV=∫_ΩψF_idV
    ∫_Ωψρu̇_idV+∫_Ω(∇ψ)ς_idV     +∫_Ωψ∂_x_ipdV=∫_ΩψF_idV.

When using integration by parts in the above, Neumann boundary condition
is imposed as 

    ς_i·dΓ=0,

where dΓ is the outward normal vector
of the tangent plane of any point on the boundary of acoustic space.
This condition means the damping force should vanish on the non-Dirichlet
boundaries. Similarly for (<ref>), multiplying
a test function vector results in 

    ∫_ΩψṗdV+∫_Ωψρ c^2(∇·u)dV=0.

To derive the ODEs, we define space discretization similar to (<ref>),
but with one more variable p, as 

    u(x,y,z,t)     =φ_1(x,y,z)·[S_1ξ(t)+S_0,1ξ_0(t)],
    
    v(x,y,z,t)     =φ_2(x,y,z)·[S_2ξ(t)+S_0,2ξ_0(t)],
    
    w(x,y,z,t)     =φ_3(x,y,z)·[S_3ξ(t)+S_0,3ξ_0(t)],
    
    p(x,y,z,t)     =φ_4(x,y,z)·[S_4ξ(t)+S_0,4ξ_0(t)].

And similar to (<ref>), we define convenience matrices
to express that 

    u_i=P_iξ+P_0,iξ_0, ∇ u_i=Q_iξ+Q_0,iξ_0, 
       P=[[ P_1; P_2; P_3 ]]=[[ φ_1^⊤S_1; φ_2^⊤S_2; φ_3^⊤S_3 ]], P_0=[[ P_0,1; P_0,2; P_0,3 ]]=[[ φ_1^⊤S_0,1; φ_2^⊤S_0,2; φ_3^⊤S_0,3 ]],
       Q=[[ Q_1; Q_2; Q_3 ]]=[[ (∇φ_1)^⊤S_1; (∇φ_2)^⊤S_2; (∇φ_3)^⊤S_3 ]], Q_0=[[ Q_0,1; Q_0,2; Q_0,3 ]]=[[ (∇φ_1)^⊤S_0,1; (∇φ_2)^⊤S_0,2; (∇φ_3)^⊤S_0,3 ]],
       P_4=φ_4^⊤S_4, P_0,4=φ_4^⊤S_0,4, Q_4=(∇φ_4)^⊤S_4 ,Q_0,4=(∇φ_4)^⊤S_0,4.

Then similar to (<ref>), we find 

    ∑_j=1^3A_ij∇ u_j=A_iQξ+A_iQ_0ξ_0,
       ∑_j=1^3B_ij∇ u_j=B_iQξ+B_iQ_0ξ_0,
       ∂_x_ip=∂_x_iφ_4·S_4ξ+∂_x_iφ_4·S_0,4ξ_0,

so that space partial derivatives can be expressed as linear transformations
of the DOF vector ξ(t) or its time derivatives.
It now sufficies to derive the system of second-order ODEs similar
to (<ref>) (<ref>).
We use P_i^⊤ as a vector of test functions
for the i th weak form equation in (<ref>)
and use P_4^⊤ to test the weak form (<ref>),
which would yield 4 groups of equations, each group having N equations.
Summing the 4 groups into 1 group, we derive a system of N coupled
first-order ODEs as 

    Mξ̇(t)+Kξ(t)=f(t),

where 

    M=∫_Ω(ρP^⊤P+P_4^⊤P_4)dV,
       M_0=∫_Ω(ρP^⊤P_0+P_4^⊤P_0,4)dV,
       K=∫_Ω[Q^⊤(μ_1A+μ_2B)Q+ρ c^2P_4^⊤HQ+P^⊤Q_4]dV,
       K_0=∫_Ω[Q^⊤(μ_1A+μ_2B)Q_0+ρ c^2P_4^⊤HQ_0+P^⊤Q_0,4]dV,
       f(t)=∫_ΩP^⊤FdV-M_0ξ̇_̇0̇(t)-K_0ξ_0(t).

Here M and K may not be as meaningful
as the mass and stiffness matrices in the second-order ODEs case.
The initial conditions are 

    u|_t=0=0, p|_t=0=0.


The acoustic space is often finite, which means boundary conditions
are essential to shape the solutions of acoustic equations. Let Γ^(c)∈ℝ^3
be the acoustic space's whole solid boundary, of which Γ_1^(c),Γ_2^(c)∈Γ^(c)
are the room boundary (grounds, walls, ceilings) and piano soundboard
parts respectively, as shown in figure <ref>. In <cit.>,
the walls were assumed rigid, which can produce sound reflection effects
but may miss sound absorption effects. This motivates our introduction
of vibratory acoustic boundaries to account for various phenomena
like reflection, scattering, absorption, diffusion, resonance etc.



 §.§ Model for the room


In this section, we delve into the specification of room boundary
models, which serve as the cornerstone for the subsequent analysis
of solid-air coupling dynamics. For the sake of conciseness, we adopt
a simplified representation of the room environment as a shoebox-shaped
3D space accommodating a piano. It is worth noting, however, that
this room model remains versatile and can be readily adapted to accommodate
irregular room geometries. The room's boundaries—encompassing floors,
walls, and ceilings—play pivotal roles in shaping the acoustic
behavior within, acting as both reflectors and absorbers of sound
waves. In the context of physical modeling, these boundaries are treated
as 3D elastic materials, akin to the piano's soundboard model. This
methodological alignment facilitates the straight application of the
general 3D elasticity model and the soundboard model to the room,
obviating the need for reformulations of the governing equations and
weak forms.

The primary specialization within this framework lies in the stipulation
of boundary conditions of room barriers. As highlighted red in figure
<ref>, the inner side of ground, wall and ceiling (blue
color), defined as Γ_1^(f), is the room material's interface
with the air; the underside of ground material (red color), defined
as Γ_2^(f), is where the Dirichlet boundary condition
applies. Firmly rooted in the earth, the ground material should experience
zero displacement along Γ_2^(f). This condition is written
as 

    u^(f)|_x∈Γ_2^(f)=0.

In case the inner concrete layers of walls and ceilings are considered
rigid, displacements on the surfaces of them should be zero too. Another
specialization is that prestress need not be considered for room barriers.
Even if it exists, the dynamic deformation may not be large enough
to induce significant effect of prestress on vibration.



§ MODEL FOR COUPLING BETWEEN DIFFERENT PARTS OF THE PIANO




 §.§ Model for hammer-string coupling




The piano hammer positioned below the string functions the excitation
of string vibration, which is known to be a highly nonlinear process
<cit.>. We consider two main parts of the hammer
relevant to this excitation process: the wooden hammer shank that
can be approximated as non-deformable; the hammer felt that is deformable,
impacting and exerting force on the string. Normally, the hammer moves
in a circle when triggered by piano player's key action, striking
the string at certain point in during a very short period of time,
and is also pushed back by the vibrating string.



  §.§.§ One hammer striking one string


For the case of one hammer striking one string, the study follows
the 0D nonlinear hammer-string interaction model in <cit.>,
with some refinements: for the case of one hammer striking one string,
the rotational movement of hammer <cit.> and
the horizontal interaction force are taken into account. 

As shown in figure <ref>, xyz is the coordinate system
of piano string, where the x axis of the string system has an angle
α_0 to the ground. To better describe the hammer shank motion
which is rotation around point P_1, we will not often the use
the xyz system. Instead, we define a static coordinate system x'y'z'
with origin P_1, y' axis the same as y axis, z' axis
with negative direction the same as gravity, and x' axis perpendicular
to the y'O'z' plane. The motion of hammer shank as a rigid body
can then be described as rotating around the y' axis and not moving
in the y' direction. We can thus use a single variable θ(t),
the angle of P_1P_3 to the x' axis, to
fuuly describe the hammer shank motion. The hammer shank drives the
overall motion of the shank head P_2 and felt head P_4.
Nevertheless, when the felt head is in contact with the string's Q_0
point, it undergoes a compression ϑ=[ϑ_1,ϑ_2,ϑ_3]^⊤
from P_4 to P_6, which also overlaps with the dynamic string
point Q_1, which also contributes to its motion. To perform a
reasonable orthogonal decomposition of this compression, we define
a dynamic coordinate system x”y”z” with origin P_3, x”
axis with positive direction P_5Q_0, y”
axis with positive direction the same as y axis, and z” axis
with positive direction P_5P_2. The hammer
compression can then be decomposed onto the x”,y”,z” axes as
marked as marked ϑ_1,ϑ_2,ϑ_3 respectively
in figure <ref>. 

Denote the displacement of string point Q_0 with coordinate x_0^(a)
in the xyz system as u_0^(a)=u^(a)(x_0^(a),t).
Denote the lengths of Q_0P_1, P_0P_1, P_1P_3,
P_2P_3, P_1P_2, P_2P_4 as L, L_0, L_1,
L_2, L_3, L_4 respectively, and the angle of P_1Q_0
to the x' axis as α_1. Our geometric analysis finds the
coordinates of Q_1 (equivalent to P_6 in case of compression)
in the x'y'z' system and x”y”z” system, denoted x_1^(a')
and x_1^(a”), are 

    x_1^(a')=r_0+R_0u_0^(a), x_1^(a”)=r+R(θ)x_1^(a'),
       r_0=[[ Lcosα_1;       0; Lsinα_1 ]], R_0=[[  cosα_0       0 -sinα_0;       0       1       0;  sinα_0       0  cosα_0 ]],
       r=[[ -L_1;    0; -L_2 ]], R(θ)=[[  cosθ     0  sinθ;     0     1     0; -sinθ     0  cosθ ]].

And converting x_1^(a”) to x_1^(a')
results in 

    x_1^(a')=s(θ)+S(θ)x_1^(a”),
       s(θ)=[[ L_1cosθ-L_2sinθ;               0; L_1sinθ+L_2cosθ ]], S(θ)=[[  cosθ     0 -sinθ;     0     1     0;  sinθ     0  cosθ ]].


As previously discussed, we assume Q_0 and P_2 is non-deformable
whereas P_4 is deformable. This means when the hammer is in contact
with string, P_4 moves to P_6; when contact is absent, the
position of P_4 is determined by the shank but not the string.
The condition to satisfy that the hammer is in contact with string
should be z_1^(a”)≤ L_4, viz. the z” direction distance
between string point and shank head is less than the static thickness
of hammer felt. The hammer felt compression is then 

    ϑ={[ x_1^(a”)-L_4, z_1^(a”)≤ L_4;             0, z_0^(a”)>L_4 ].

where L_4=[0,0,L_4]^⊤; negative value of
ϑ_i means the hammer felt is compressed towards the negative
direction of x_i” axis, and vice versa. Due to compression,
the hammer felt exerts an interaction force F^(a”)=[F_1^(a”),F_2^(a”),F_3^(a”)]^⊤
(in x”y”z” system) on the string. Reciprocally, the hammer felt
suffers -F^(a”) from the string according to Newton's
third law[Note that we ignore here the string's normal and tangent stress and
prestress on the interaction surface exerting on the hammer felt (as
well as the relevant strain on the string side), which should already
be zero per the Neumann boundary condition specified in (<ref>)
if DOFs are given to string displacements on this surface. Even if
they should not be zero, it is relatively acceptable to ignore them
as they are probably less contributive than the hammer's compression
force. Treating them as non-zero would require imposing Dirichlet
boundary conditions on some relevant DOFs of the string's side, which
is a computational challenge.]. Based on <cit.>, the hammer-string interaction
force in the x”_i axis is modelled as a nonlinear function of
compression as 

    F_i^(a”)=-sgn(ϑ_i)[k_i|ϑ_i|^p_i+r_ik_i∂_t(|ϑ_i|^p_i)], i=1,2,3,

where k_i is the stiffness of hammer felt; p_i is a positive
exponent accounting for nonlinearity; r_i is the relaxation coefficient
accounting for the hammer's hysteretic and dissipative behaviour;
function sgn(·) returns the sign of a real number.
Converting the interaction force to x'y'z' and x'y'z' coordinates
yields F^(a')=S(θ)F^(a”)
and F^(a)=R_0^⊤F^(a').

To couple the hammer force with string motion, we need equations in
the x'y'z' system governing the motion of hammer shank, a rigid
body. Denote the coordinates of shank mass centre P_0, shank
rotation centre P_1, shank head P_2 in the x'y'z' system
as x_0^(d'), x_1^(d'), x_2^(d')
respectively. It can be found that x_0^(d')=[L_0cos(β_0+θ),0,L_0sin(β_0+θ)]^⊤
and x_2^(d')=s(θ), where β_0
is the angle of P_1P_0 to the P_1P_3.
Define the total mass and homogenous line density of the hammer shank
as m, ρ. It can be found that m=ρ(L_1+L_2), and
the axis-free position of center of mass P_0 projected on P_1P_3
and P_2P_3 are (1/2L_1^2+L_1L_2)/(L_1+L_2)
and 1/2L_2^2/(L_1+L_2) respectively. It is obvious
that the shank as a rigid body suffers these forces: string reaction
force F^(d')=-F^(a') at P_2;
a rotation constraint force at P_1, which has zero torque with
respect to the y' axis; gravity [0,0,-mg]^⊤ at P_0,
where g is the scalar gravitational acceleration. It follows that
the total torque with respect to the y' axis is mgL_0cos(β_0+θ)+s_2(θ)·F^(d')
(only 1 dimension of torque is needed here), where 

    s_2(θ)=[[  L_1sinθ+L_2cosθ;                0; -L_1cosθ+L_2sinθ ]];

the moment of inertia with respect to the y” axis, denoted I,
is 

    I=∫_0^L_1ρ l^2dl+∫_0^L_2ρ(l^2+L_1^2)dl=ρ(1/3L_1^3+1/3L_2^3+L_1^2L_2).

Applying Newton's second law for rotation, the differential equation
of hammer shank motion is 

    Iθ̈=-μθ̇+mgL_0cos(β_0+θ)+s_2(θ)·F^(d'),

where -μθ̇ is the simplified damping force. Here t=0
is the moment when the felt head first contacts the string and compression
is still zero, and initial angle and angular velocity at this time
should be known. 

From the aforementioned deductions, we can abstract the hammer-string
interaction force F^(e) as a nonlinear operator
ℱ(u_0^(a),θ):(ℝ^3,ℝ)→ℝ^3.
This force should be added to the rhs of (<ref>) as a non-conservative
force. Then, (<ref>) and (<ref>) forms
2 sets of coupled equations for 2 sets of variables (u^(a),θ),
which provides a foundation for obtaining a solution theoretically.
Due to the highly nonlinear nature of hammer-string coupling, particularly
in that Taylor series approximation may not work well for potentially
non-smooth functions in (<ref>), it would be inappropriate
to use the common perturbation method. In section <ref>
we shall introduce an explicit time discretization method to efficiently
solve the hammer-string coupling.



  §.§.§ One hammer striking two or three strings


The case of one hammer striking multiple, say three strings can simply
be treated it as if three independent felt heads were striking their
corresponding strings. Three independent interaction forces F_i=1,2,3^(e')
are computed from three independent compressions. The string force
exerted on the shank is the sum F^(d')=-∑_i=1^3F_i^(e'),
assuming the 3 compression forces apply to the same point of shank
head. However, this approach may lose the dynamic interaction between
different striking points of the hammer felt.

A 3D hammer model may be more accurate in capturing the interplay
of different hammer striking points and the different positions of
compression forces. The hammer felt is now modelled as a 3D elastic
material with space discretization. Two kinds of spacial boundaries
consist in the hammer felt: one contains the felt's 3 potential contact
points with the 3 strings, represented as P_4,i with coordinates
x_4,i=1,2,3^(e”); another is the contact surface
Γ_1^(e”)∈ℝ^2 with the hammer shank. The
hammer shank is considered as a rigid 3D object, which means no inner
space discretization is needed for it. Above all, one needs to pay
attention to the choice of coordinate system for the 3D geometry of
felt. A recommended choice here is the dynamic x”y”z” system,
which eliminates the felt's rigid body movement component and is thus
suitable for FEM computation. Also, displacement fields u^(e”)(x^(e”),t)
for hammer felt and u^(d”)(x^(d”),t)
for hammer shank need to be defined.

Derivation of the felt compression dynamics is now based on the 3D
elastic material model and soundboard model previously introduced,
without the need for modelling prestress. We first consider the boundary
conditions needed for the felt model. Denote the coordinate of Q_1,i
in the x”y”z” system as x_1,i^(a”), which
can be computed similar to (<ref>).
To represent hammer felt compression, Dirichlet boundary condition
should be imposed on these 3 contact points as 

    u^(e”)(x_4,i^(e”),t)={[ x_1,i^(a”)-x_4,i^(e”), z_1,i^(a”)≤ L_4;                      0, z_1,i^(a”)>L_4 ].,

where z_1,i^(a”)≤ L_4 means contact is present and otherwise
not. This determined displacement would occur on the rhs of system
ODEs like (<ref>) as a source term. This means
without contact, not only P_4,i but also the entire felt will
not experience deformation; with contact, the position of P_4,i,
i.e. P_6,i, is the same as Q_1,i, triggering motion and
compression of the hammer felt. As for another boundary Γ_1^(e”)
which interfaces the shank, Dirichlet boundary condition should be
imposed on as 

    u^(e”)(x^(e”),t)=0, x^(e”)∈Γ_1^(e”),

because the rigid shank should have no deformation.

We then consider the dynamics of string and shank impacted by the
felt. At the contact point with string, the felt's nonzero surface
forces G^(e”)(x_4,i^(e”),t) should
be transmitted to the string as 

    F_i^(a)=-[[    0 Z_12 Z_13; Z_12 Z_22 Z_23; Z_13 Z_23 Z_33 ]][[ 1; 1; 1 ]],
       { Z_ij} =R_0^⊤R(θ)^⊤G^(e”)(x_4,i^(e”),t)R(θ)R_0,

where Z_11 is not transmitted considering that the hammer felt
seems to have no direct contact with the string's surface whose normal
is in the x axis (longitudinal direction). As for the felt's surface
force exerting on the shank along the boundary Γ_1^(e”)=Γ_1^(d”),
previous practice of directly transmitting the felt-string interaction
force to the shank is not applicable here due to the 3D nature of
felt and shank. For a point x^(e”)=x^(d”)
in Γ_1^(e”), with n^(e”)(x^(e”))
as the outward normal of its tangent plane, the felt would exert a
surface force 

    F^(d”)(x^(d”),t)=-G^(e”)(x^(e”),t)n^(e”)(x^(e”))

 on the shank. Then the total contribution of felt force to the torque
of shank with respect to the y” axis is 

    T_1=∫_Γ_1^(d')(F_1^(d')z-F_3^(d')x)dV

where F_i^(d') is the i th entry of vector F^(d')(x^(d'),t)=S(θ)F^(d”)(R(θ)x^(d'),t),
and Γ_1^(d') is converted from Γ_1^(d”) using
(<ref>). The total contribution
of gravity to the torque of shank with respect to the y” axis
is T_0=mgx_0^(d') where x_0^(d') is the
coordinate of P_0 depending on θ. The moment of inertia
with respect to the y” axis is 

    I=∫_Ω^(d')ρ(x^2+z^2)dV,

where Ω^(d')∈ℝ^3 is the shank's volume domain.
Applying Newton's second law for rotation, the differential equation
of hammer shank motion is 

    Iθ̈=-μθ̇+mgx_0^(d')+∫_Γ_1^(d')(F_1^(d')z-F_3^(d')x)dV.


Now the whole model of 3D hammer-string coupling has been established.
The felt-string interaction force F_i^(e) in (<ref>)
is computed from the 3D hammer felt model with Dirichlet boundary
conditions (<ref>) (<ref>),
and added to the rhs of (<ref>) as a non-conservative force.
Then, (<ref>) and (<ref>) forms 2 sets
of coupled equations for 2 sets of variables (u^(a),θ),
which provides a foundation for obtaining a solution theoretically.
Despite the linearity of 3D elasticity model, the felt's force on
the string is still nonlinear because of the interaction judging condition.



 §.§ Model for string-soundboard coupling




Piano strings are coupled to the soundboard's bridge part, which terminates
and transmits string vibration to the soundboard. Figure <ref>
visualizes this coupling at the bridge from 3 perspectives, where
xyz and x'y'z' are the coordinate systems for soundboard and
string respectively. Converting a vector (not a point) v
in the xyz system into v' in the x'y'z' system,
according to the angles α and β marked in figure <ref>,
yields v'=Rv where R=R_1R_2
and 

    R_1=[[  cosα     0  sinα;     0     1     0; -sinα     0  cosα ]], R_2=[[  cosβ  sinβ     0; -sinβ  cosβ     0;     0     0     1 ]]

are the rotation matrices. The orthogonal property of rotation matrix
yields v=R^⊤v'.

We conjecture from observations that the string-soundboard coupling
is achieved via several mechanisms. The first mechanism is “bridge
hump coupling”, as illustrated in the side perspective of figure
<ref>. The middle part of bridge is constructed a bit
higher than the agraffe, the starting point of the speaking string.
As a result, the string experiences some upward pressure from the
bridge and the bridge experiences some downward pressure from the
string, both statically and dynamically. The second mechanism is “bridge
pin horizontal coupling”, as illustrated in the top perspective
of figure <ref>. Two bridge pins for one string are drilled
into the bridge at positions that would bend the string a bit horizontally,
bringing mutual horizontal tension between the string and bridge,
and restricting the string's horizontal movement. Made of hard metal
like steel, brass or even titanium, the bridge pins can be treated
as rigid bodies with zero strain and stress. The third mechanism is
“bridge pin vertical coupling”, as illustrated in the front perspective
of figure <ref>. The two bridge pins lean towards their
respective string sides, blocking the string from moving upwards beyond
the pins. Also, the notches of bridge ensure that bridge hump coupling
and bridge pin coupling occur at almost the same position, viz. the
front bridge pin (the one closer to the agraffe). 

On the string's side, the coupling point is defined as x_1^(a)=[L_1^(a),0,-r^(a)]^⊤,
the lowest point of string at the front pin. While other relevant
points for coupling may be [L_1^(a),-r^(a),0]^⊤, the
closest point of string to the bridge pin, the string's displacements
and surface forces at these points are identical, assuming u^(a)
does not vary with y^(a) or z^(a) at x^(a)=L_1^(a)
(no rotation) for simplicity[This is also true for the string's central point x_1^(a)=[L_1^(a),0,0]^⊤
previously defined.]. The coordinate of this coupling point on the soundboard's side is
denoted x_1^(b).



  §.§.§ Static coupling


For coupling before motion is initiated, the string and soundboard
should have non-zero prestress fields that attain static balance.
It is sufficient to specify that the string prestress only acts perpendicular
to its cross-sections and is constant through its speaking length.
This means in the string's tension matrix T^(a),
T_11 is a constant and other entries are zero, immediately satisfying
the static balance. For the soundboard, a full space dependent tension
matrix T^(b) with 6 unique elements is necessary.
The static balance condition can then be expressed as 3 partial differential
equations 

    ∇·T^(b)=0,

which can be solved numerically by FEM using nodal surface forces
as time-independent DOFs. Nonetheless, continuity between string and
soundboard tension fields is required. To write the continuity equations,
we define an intermediate coordinate system x”y”z” resulting
from rotating the x'y'z' system by R_1^⊤
or rotating the xyz system by R_2, as marked
blue in figure <ref>. At the coupling point, the string's
tension on surfaces with outward normals in the y” and z”
axes exert on the bridge, whereas that in the x” axis exerts on
the farther hitch pins that seem to not connect the soundboard. This
leads to one Dirichlet boundary condition for each string coupled
to the soundboard as 

    (R_2T^(b)(x_1^(b))R_2^⊤)_[22,33,12,13,23]=(R_1^⊤T^(a)(x_1^(a))R_1)_[22,33,12,13,23],

where subscript _[i_1,i_2,...] means taking vector or matrix
elements at indices i_1,i_2,... to form a vector. This condition
can be deduced into 5 equations, leaving 1 DOF for the coupling point. 



  §.§.§ Dynamic coupling


For coupling in motion, the string transmits its surface forces (dynamic
prestress, stress and damping force) to the soundboard. At coupling
point, the string's total surface force, excluding the static prestress
which is already accounted for on the soundboard's side, is G^(a)(x_1,t).
Similar to static coupling, only forces on surfaces with outward normals
in the y” and z” axes are transmitted to the soundboard at
the coupling point. This leads to 

    F^(b)(x^(b),t)=δ(x^(b)-x_1^(b))(R_2^⊤[[    0 Z_12 Z_13; Z_12 Z_22 Z_23; Z_13 Z_23 Z_33 ]]R_2)[[ 1; 1; 1 ]],
       { Z_ij} =R_1^⊤G^(a)(x_1^(a),t)R_1,

where δ(x) is a 3D Dirac delta function indicating
a point load. The above equations treats the string's surface force
as a non-conservative force input to the soundboard system.

Another coupling condition arising from observation is that the string
and soundboard should have the same y” and z” direction displacements[Here displacement continuity should be equivalent to velocity continuity,
since both the string's and soundboard's displacements are based on
a [0,0,0]^⊤ point.] at the coupling point, written as 

    (R_1^⊤∂_tu^(a)(x_1^(a),t))_[2,3]=(R_2∂_tu^(b)(x_1^(b),t))_[2,3].

The reason for this is that the the aforementioned three coupling
mechanisms constitute support for the string's surfaces with outward
normals in the y” and z” axes at the coupling point. This
may not hold for the x” axis, but if it is needed[Maybe in case of high enough viscosity or friction, the string can
not slip more than the soundboard in the x” direction. Restricting
this movement may also make the solution of string's PDE more stable.
This needs to be tested.], the subscript _[2,3] in above equations can be removed to impose
a stronger displacement continuity condition. Anyway, restriction
of the string's x” axis (longitudinal) motion is always present
at the farther hitch point which may not connect the soundboard, see
(<ref>). The same-displacement condition also
provides a basis for establishing the string's displacement at the
coupling point (Dirichlet boundary condition), so that the string's
PDE (<ref>) has an appropriate solution.



  §.§.§ More discussions on coupling


String-soundboard coupling exhibits complex mechanisms that we find
challenging to discover and describe. Some of these mechanisms that
we observe but not covered in the aforementioned coupling model are
discussed in this subsection. These discussions may not cover much
details of computation due to their complexities.

Firstly, we only specified the coupling point at the bridge's front
pin but ignored the rear pin. As coupling at the rear pin position
potentially exists, the string in our model may need to be extented
in length to cover the segment between the front pin position x_1=[L_1^(a),0,-r^(a)]^⊤
and the rear pin position x_2=[L_2^(a),0,-r^(a)]^⊤
(we call it “this segment” is this paragraph). At this segment,
the horizontal bending of string may be accouted for in the prestress,
but not necessarily represented in the geometrical volume for simplicity.
In static state, the prestress field coupling should be computed for
this segment paying attention that only at the two pin points are
y” axis surface forces coupled. In dynamic state, transmission
of the string's surface forces should be computed for this segment,
paying attention that only at the two pin points are y” axis surface
forces transmitted. A maybe more accurate computation of surface force
transmission is, that at the two pins only when the string's y”
axis surface forces are towards the front or rear pins should they
be transmitted, and that at between (excluding) the two pins only
when the string's z” axis surface forces are upwards should they
be transmitted. This may lead to a highly nonlinear function akin
to the case of hammer-string interaction (<ref>).
For displacement coupling at this segment, at the two pins when the
string's y” direction displacements should not exceed the pins,
whereas the z” direction displacements are fully coupled to the
pins; at between (excluding) the string's z” direction displacements
should not be under the soundboard plane, whereas the y” direction
displacements are unrestricted.

Secondly, we considered coupling as occurring at a point, but it may
actually occur in a small contact surface (we call it “this surface”
is this paragraph). In such a case, the Dirichlet boundary condition
of prestress field coupling should be specified over this surface;
the dynamic surface force transmission should be treated as surface
load rather than point load; the Dirichlet boundary condition of displacement
coupling should also be specified over this surface.

Thirdly, we still lack geometric details regarding the string's notable
vertical and horizontal bending at the bridge. This bending changes
the string's longitudinal direction, and thereby the 3 directions
of vibration. Though string vibration should be terminated at the
bridge pin, it seems only at the hitch pin that longitudinal vibration
is fully restricted. Therefore, the string's segment from front bridge
pin to hitch pin may require investigation, which may concern the
duplex scale phenomena <cit.>. We can design
a multi-segment geometric model for the string, each segment having
differently rotated constitutive matrices. Also, the position-dependent
static tension can be specified parallel to the central line of each
string segment.

Finally, we ignored the soundboard's normal and tangent surface forces
exerting on the string at coupling point. This is similar to the case
of hammer-string coupling where the string's surface forces exerting
on the hammer are ignored. This is mainly for practical considerations,
as we would want to solve the soundboard's equations for only once.
Since we did not impose Dirichlet boundary condition for the soundboard
at each coupling point with the around 200 strings, Neumann boundary
conditions apply here restricting the soundboard's conservative surface
forces to be zero at coupling points. It may be acceptable given the
intuitive feeling that “string → soundboard” transmission should
dominate “soundboard → string” transmission. If the latter is
essential, some DOFs may be need to be removed from the coupling point
on the soundboard's side, which is a challenge for imposing Dirichlet
boundary condition reasonably for both the string and the soundboard.



 §.§ Model for soundboard-air and room-air coupling


Modelling soundboard-air and room-air coupling is crucial for arriving
at the final digital audio sound to the listeners. Generally in the
context of solid-fluid coupling dynamics, the interaction surface
should have continuous normal and tangential velocities, as well as
continuous normal and tangential surface forces on both sides <cit.>. 

As per section <ref>, room-air
coupling occurs on surface Γ_1^(f) (room side) and Γ_1^(c)
(air side). Define coordinate transformation x^(f)=r_1+R_1x^(c)
from air coordinates to soundboard coordinates, where r_1
and R_1 are the shift vector and rotation matrix.
We impose Dirichlet boundary condition of velocity continuity that

    u^(c)(x^(c),t)=R_1^⊤u̇^(f)(x^(f),t), x^(c)∈Γ_1^(c), x^(f)=r_1+R_1x^(c).

We also specify the contribution of air surface force to the room
equation's source term on the rhs of (<ref>) as 

    F_1^(f)(x^(f),t)=-G^(c)(x^(c),t)n^(c)(x^(c)), x^(f)∈Γ_1^(f), x^(c)=-R_1^⊤r_1+R_1^⊤x^(f),

where n^(c)(x^(c)) is the outward
normal of the tangent plane of x^(c).

As per sections <ref> and <ref>,
soundboard-air coupling occurs on surface Γ_2^(b) (soundboard
side) and Γ_2^(c) (air side). Define coordinate transformation
x^(b)=r_2+R_2x^(c)
from air coordinates to soundboard coordinates, where r_2
and R_2 are the shift vector and rotation matrix.
We impose Dirichlet boundary condition of velocity continuity that

    u^(c)(x^(c),t)=R_2^⊤u̇^(b)(x^(b),t), x^(c)∈Γ_2^(c), x^(b)=r_2+R_2x^(c).

We also specify the contribution of air surface force to the soundboard
equation's source term on the rhs of (<ref>) as 

    F_1^(b)(x^(b),t)=-G^(c)(x^(c),t)n^(c)(x^(c)), x^(b)∈Γ_2^(b), x^(c)=-R_2^⊤r_2+R_2^⊤x^(b),

where n^(c)(x^(c)) is the outward
normal of the tangent plane of x^(c).



§ NUMERIC SCHEMES




 §.§ Modal superposition method for solving coupled ODEs




  §.§.§ Modal transformation of second-order ODEs


For solving coupled second-order ODEs like (<ref>),
time discretization using finite-difference is a direct approach.
However, given that damping matrix is diagonalizable by the mass and
stiffness matrices, decoupling and dimension reduction of the system
by means of modal superposition is preferred. To do so, we first define
the following eigenvalue problem 

    Kϕ_i=λ_iMϕ_i⟺M^-1Kϕ_i=λ_iϕ_i (i=1,...,N),

and the eigen decomposition 

    K=MΦΛΦ^-1⟺M^-1K=ΦΛΦ^-1,

where λ_i is a complex eigenvalue, ϕ_i
is a N×1 complex eigenvector, Λ=diag(λ_1,...,λ_N)
is a diagonal matrix of eigenvalues, Φ=[ϕ_1,...,ϕ_N]
is a N× N matrix of eigenvectors. Note that since both M
and K are sparse matrices with large dimensions, it
is preferrable to solve the generalized eigenvalue problem on the
left side of (<ref>), avoiding
explicit inversion of M that would otherwise result
in a large dense matrix. We can leverage existing softwares of sparse
eigensolvers like Arpack and FEAST to solve the generalized eigenvalue
problem. Nevertheless, for solving ODEs of the acoustic system with
large DOFs, it may be unwise to store all rows of eigenvectors at
the same time which may otherwise lead to memory overload. A re-implementation
of existing sparse eigensolver algorithms may be desired so as to
store partial rows or columns of eigenvectors in a “rolling” way.

Since eigenvectors are linearly independent, that is to say Φ
is invertible, we can write the solution in the form ξ(t)=Φq(t)
where q(t) is an N×1 vector we call as modal
DOFs. Note here ξ(t) should be real (in the complex
domain), but q(t) is complex because Φ^-1
is complex. Then (<ref>) can be decoupled as 


    MΦq̈(t)+2μKΦq̇(t)+KΦq(t)     =f(t)
    Φ^HMΦq̈(t)+2μΦ^HKΦq̇(t)+Φ^HKΦq(t)     =Φ^Hf(t)
    M̅q̈(t)+2μM̅Λq̇(t)+M̅Λq(t)     =Φ^Hf(t)
    q̈(t)+2μΛq̇(t)+Λq(t)     =p(t)


 where 

    M̅=Φ^HMΦ, p(t)=M̅^-1Φ^Hf(t);

are called the modal mass matrix and the modal force vector respectively.
Through eigen decomposition, the N coupled ODEs have now been transformed
into N uncoupled ODEs. 

It is well-understood that the analytical solution of the i th
uncoupled second-order ODE without the rhs nonhomogenous term is a
sinusoidal signal with a single eigenfrequency positively correlated
to the magnitude of i th eigenvalue. Assuming eigenvalues in Λ
are ascendingly sorted by their real magitudes, we can select only
the lowest M eigenvalues and discard the rest, because most human
ears are insensitive to eigenfrequencies above a certain threshold
(often 10 kHz). This is also for practical considerations that the
number of DOFs is often too large for numeric computation, making
dimension reduction desirable. Consequently, we obtain a M×1
vector q'(t) with only the first M entries of q(t),
and the M+1 to N entries are all approximated as zeros and discarded.
Subsituting q(t) by q'(t) (more exactly,
q'(t) should be padded N-M zeros at the end) in
(<ref>), and omitting the last N-M
equations and unknowns, yields 


    S'q̈'(t)+2μS'Λ'q̇'(t)+S'Λ'q'(t)     =Φ'^⊤f(t)
    q̈'(t)+2μΛ'q̇'(t)+Λ'q'(t)     =p'(t)


 where Λ' has dimension M× M containing
the first M eigenvalues and Φ' has dimension
N× M containing the first M eigenvectors. The large numbe
of DOFs is now approximated as the linear combination of a smaller
number of modal DOFs as ξ(t)≈Φ'q'(t).
The reduced modal mass and modal force are 

    M̅'=Φ'^HMΦ', p(t)=M̅'^-1Φ'^Hf(t).

For around 2500 modes of the soundboard system as estimated in <cit.>,
the fully dense modal mass matrix (128 bit complex type) would require
about 95 MB memory which is normally acceptable in both storage and
computation aspects. If the mass and stiffness matrices are both symmetric,
then eigenvectors can be normalized so that the modal mass matrix
equals to an identity matrix, making its inversion easier. However,
this benefit is not achievable in our piano model because the constitutive
matrix would lose its symmetry in the presence of prestress.

Now we consider the analytical solution of the i th uncoupled nonhomogenous
ODE 

    q̈(t)+2μλq̇(t)+λ q(t)=p(t),

where the subscripts _i and superscripts ' are dropped for
convenience. Applying Laplace transform to this equation yields 

    s^2Q(s)+2μλ sQ(s)+λ Q(s)     =P(s)+Q_0(s),
    
    s^2G(s)+2μλ sG(s)+λ G(s)-G_0(s)     =1
    g̈(t)+2μλġ(t)+λ g(t)     =δ(t)

where s is a complex variable, δ(t) is the Dirac delta
function; the initial parts of Laplace transforming derivatives are
defined as 

    Q_0(s)     =sq(0)+q̇(0)+2μλ q(0)
    
    G_0(s)     =sg(0)+ġ(0)+2μλ g(0);

the well-known Green's function (frequency domain), the response to
a unit impulse, is defined as 

    G(s)=Q(s)[1+G_0(s)]/P(s)+Q_0(s)⇒ Q(s)=G(s)[P(s)+Q_0(s)]/1+G_0(s).

To solve g(t), we first specify that response should not exist
during zero and negative time, viz. g(t)=ġ(t)=0 for t≤0;
then notice when t>0, (<ref>) becomes
homogenous with solution 

    g(t)=C_1exp(z_1t)+C_2exp(z_2t), z_1,z_2=-μλ±√((μ^2λ-1)λ).

To determine complex constants C_1, C_2, we first notice
that g(t) should be continuous at t=0 in the presence of second
derivative in (<ref>), thus 

    g(0^+)=g(0)⇒ C_1+C_2=0.

Then, integrating (<ref>) over (-∞,+∞)
yields 

    ∫_-∞^+∞[g̈(t)+2μλġ(t)+λ g(t)]dt=∫_-∞^+∞δ(t)dt
    ⇒   ∫_0^-^0^+g̈(t)dt=ġ(0^+)-ġ(0^-)=ġ(0^+)=1
    ⇒    C_1z_1+C_2z_2=1,

which is the so-called jump discontinuity condition for first derivative;
the continuity condition for the zero-order has been applied again
here. From (<ref>) (<ref>)
the complex constants are found to be 

    C_1=1/z_1-z_2, C_2=1/z_2-z_1.

Subsituting g(0)=0 and ġ(0)=0 into (<ref>)
yields G_0(s)=0. Then from (<ref>)
we have 

    Q(s)     =P(s)G(s)+[q̇(0)+2μλ q(0)]G(s)+sq(0)G(s)
    
    q(t)     =p(t)*g(t)+[q̇(0)+2μλ q(0)]g(t)+q(0)ġ(t),

where * is the convolution operator over [0,+∞). For the
general case of initial conditions q(0)=q̇(0)=0, the solution
reduces to q(t)=p(t)*g(t). It is now clear that the nonhomogenous
ODEs can be solved by convolving the modal force p(t) (source signal)
with the Green's function g(t) (response signal), which can be
efficiently computed using fast Fourier transform (FFT) convolution
in the frequency domain. However, the analytical solution may not
be actually useful or efficient in the presence of coupling between
systems and the explicit time discretization scheme introduced in
section <ref> will be an alternative.
Nevertheless, it provides an understanding of the characteristics
of vibration modes.



  §.§.§ Modal transformation of first-order ODEs


Solving coupled first-order ODEs like (<ref>) by means
of modal superposition is similar to the second-order case, and we
only discuss some particularities here. Firstly, the eigenvalue problem
is defined in the same way for M (now for first-order)
and K. The decoupled ODEs write 


    M̅q̇(t)+M̅Λq(t)     =Φ^Hf(t)
    q̇(t)+Λq(t)     =p(t)


 where 

    M̅=Φ^HMΦ, p(t)=M̅^-1Φ^Hf(t)

are the modal mass and the modal force. Specific attension should
be paid to the analytical solution of the i th uncoupled nonhomogenous
ODE 

    q̇(t)+λ q(t)=p(t).

Laplace transform yields 

    sQ(s)+λ Q(s)     =P(s)+q(0),
    
    sG(s)+λ G(s)-g(0)     =1
    ġ(t)+λ g(t)     =δ(t)

where the Green's function is defined as 

    G(s)=Q(s)[1+g(0)]/P(s)+q(0)⇒ Q(s)=G(s)[P(s)+q(0)]/1+g(0).

Note that g(t)=ġ(t)=0 for t≤0 is still required, but
continuity g(0)=g(0^+) is unnecessary because the highest order
of derivative in (<ref>) is only one.
Nevertheless, continuity of ∫ g(t)dt at t=0 can
be easily satisfied because the primitive function can have an arbitrary
constant added. Then integrating (<ref>),
we can find g(0^+)=1 and the solution of Green's function 

    g(t)=-1/λexp(-λ t).

It follows that 

    Q(s)     =P(s)G(s)+q(0)G(s)
    
    q(t)     =p(t)*g(t)+q(0)g(t),

is the solution of the uncoupled nonhomogenous first-order ODE.



 §.§ Explicit time discretization for coupling between systems


Having discussed the weak form ODEs for each subsystem of the physical
piano, our question is then how a numerical treatment of coupling
between system may be achieved that attains a sensible tradeoff between
feasibility and accuracy. For most subsystems presented in the previous
sections, there exists an external force term on the rhs of PDEs.
This, along with the predetermined displacements on the boundary,
form the main contributions to the rhs source term of system ODEs. 

Nevertheless, to say “ODEs” here is actually indefensible, because
in many coupling cases discussed before the rhs source term f(t)
depends linearly or nonlinearly on the lhs unknown DOFs ξ(t)
too. In such cases, eigen decomposition of “pseudo ODEs” can only
obtain decoupled lhs operators on q(t) on the lhs,
but rhs operators on q(t) often remain coupled, not
only because of the nonlinearity of operators, but also because eigen
decomposition applies to the local subsystem but not the global system
where rhs sources come from. Hammer-string coupling is a typical example:
the string displacement depends on the felt's force, the felt's force
depends on its compression, but this compression depends on string
displacement. Another example is in soundboard-string coupling, the
soundboard displacement at the coupling point depends on the string's
surface force, the string's surface force depends on its displacement,
computing its displacement requires knowing the coupling point displacement
(Dirichlet boundary), but this displacement depends on the soundboard
displacement at the coupling point to satisfy the continuity condition. 

An explanation for these seemingly odd relations is that a theoretically
sound way would seem to do computations as if all subsystems were
a single system, rather than to separate systems and reintroduce coupling
between them. This means all DOFs are integrated into a single vector
of DOFs, all mass (or stiffness) matrices are assembled into a bigger
one, displacements at coupling positions are expressed by the same
unique DOFs behind, in order for full coupling solutions. However,
this is often unacceptable not only because of the high costs of storage
and computation, but also because the inherent heterogeneities of
different subsystems, particularly the nonlinearity of hammer felt
compression, the first-order characteristic of acoustic system, the
different damping mechanisms, may actually not cohere well and may
even be tough to control in a single system. Therefore, we shall still
adopt the framework of separate systems with mutual coupling, and
seek for numeric schemes capable of achieving a level of accuracy
as close to full coupling schemes as possible.

The time domain scheme we shall introduce here is inspired by the
idea of velocity Verlet algorithms, but customized for our case of
modal-transformed first and second order ODEs, with second-order accuracy[Alternatively, one may consider transforming time domain into frequency
domain to explore coupling from the perspective of mobility <cit.>.]. Despite being lower-order compared to the complex time schemes in
<cit.>, our scheme
may be more efficient particularly in eliminating the need to invert
or solve big or dense matrices at each step, even when the mass matrix
is non-diagonal. Due to the existence of coupling, the analytical
solutions of first and second order ODEs discussed in section <ref>
may be not applicable here. Nevertheless, system decoupling and dimension
reduction via modal decomposition will be shown useful in improving
the efficiency of time-stepping algorithms. 



  §.§.§ Time stepping of second-order ODEs


In the previously derived second-order decoupled ODEs (<ref>),
the rhs source term can rewritten as 

    q̈(t)+2μΛq̇(t)+Λq(t)=p(t)=r[t,q(t),q̇(t)],

where r incoporate all DOF-independent and DOF-dependent
contributions to the source term; DOF-dependent contributions consist
of two major sources: non-zero Dirichlet boundary conditions; non-conservative
force (but excluding damping force already on the lhs). 

To perform time discretization, we define the discrete time interval
as h=Δ t, which can be chosen as 1/44100 seconds for producing
44.1 kHz digital audio; define the n th discrete point in time
as t_n=nh. The notion of discrete points in time form the basis
of time stepping algorithms, where integrations are performed between
time steps. An integral over a small interval [x_0,x_0+Δ x]
can be approximated as area of trapeziud 

    ∫_x_0^x_0+hf(x)dx≈h/2[f(x_0)+f(x_0+h)],

which can achieve a relatively high accuracy, though f(x_0+h)
may not be known beforehand in some contexts. If we seek for a first-order
approximation of f(x_0+h), then 

    ∫_x_0^x_0+hf(x)dx≈ hf(x_0)+h^2/2f'(x_0)

is a less accurate approximation. Utilizing these integration strategies,
we integrate (<ref>) over [t_n,t_n+1] to
get 

    q̇(t_n+1)-q̇(t_n)     =-2μΛ∫_t_n^t_n+1q̇(t)dt-Λ∫_t_n^t_n+1q(t)dt+∫_t_n^t_n+1p(t)dt
    q̇(t_n+1)-q̇(t_n)    ≈-hμΛ[q̇(t_n)+q̇(t_n+1)]-h/2Λ[q(t_n)+q(t_n+1)]+hp(t_n)+h^2/2ṗ(t_n).

Also it is obvious that 

    q(t_n+1)=q(t_n)+∫_t_n^t_n+1q̇(t)dt≈q(t_n)+h/2[q̇(t_n)+q̇(t_n+1)].

The solution of (<ref>) (<ref>)
is 

    q̇(t_n+1)    ≈Z_1^-1[Z_0q̇(t_n)-hΛq(t_n)+hp(t_n)+h^2/2ṗ(t_n)],

where 

    Z_0=I-(hμ+h^2/4)Λ, Z_1=I+(hμ+h^2/4)Λ

are diagonal matrices easy to invert. (<ref>)
(<ref>) reveal that q(t_n+1),
q̇(t_n+1) can be computed with second-order
accuracy using q(t_n), q̇(t_n),
p(t_n), ṗ(t_n), but without
using unknowns of the t_n+1 step. This means our time discretization
scheme is explicit, though it remains to discover its energy identity.
The most significant accuracy loss of this scheme seems to be the
use of (<ref>) rather than (<ref>)
in approximating the integral of p(t), where would
expect the first-order approximation of p(t_n+1)
to be acceptable with small h. This actually implies that we rely
on some historical (t_n) displacement and velocity values to
predict the current (t_n+1) input sources, then the current input
sources are updated by the computed current displacement and velocity
values. This process can be repeated several times for one time step
if higher accuracy is desired, able to use (<ref>)
rather than (<ref>) in approximating the integral
of p(t).



  §.§.§ Time stepping of first-order ODEs


Recall the previously derived second-order decoupled ODEs (<ref>),
the rhs source term can rewritten as 

    q̇(t)+Λq(t)=p(t)=r[t,q(t)],

where r incoporate all DOF-independent and DOF-dependent
contributions to the source term; DOF-dependent contributions consist
of two major sources: non-zero Dirichlet boundary conditions; non-conservative
force (but excluding damping force already on the lhs). Integrating
(<ref>) over [t_n,t_n+1], the approximate
solution of q(t_n+1) can be found as 

    q(t_n+1)-q(t_n)     =-Λ∫_t_n^t_n+1q(t)dt+∫_t_n^t_n+1p(t)dt
    q(t_n+1)-q(t_n)    ≈-h/2Λ[q(t_n)+q(t_n+1)]+hp(t_n)+h^2/2ṗ(t_n)
    q(t_n+1)    ≈(I+h/2Λ)^-1[(I-h/2Λ)q(t_n)+hp(t_n)+h^2/2ṗ(t_n)].

Therefore, q(t_n+1) can be computed using q(t_n),
p(t_n), ṗ(t_n), but without
using unknowns of the t_n+1 step.



  §.§.§ Time stepping of hammer shank rotation ODE


The equation of hammer shank motion (<ref>)
can be written in a more abstract way as 

    Iθ̈=-μθ̇+T(θ),

where T(θ) incoporate the contributions of gravity and felt
force to the total torque of shank, dependent on the zero-order value
of θ. Similar to second-order system ODEs, the time discretization
of (<ref>) is found to be 

    I[θ̇(t_n+1)-θ̇(t_n)]     =-μ∫_t_n^t_n+1θ̇dt+∫_t_n^t_n+1T(θ)dt
    
    I[θ̇(t_n+1)-θ̇(t_n)]    ≈-h/2μ[θ̇(t_n)+θ̇(t_n+1)]+hT(t_n)+h^2/2Ṫ(t_n)
    θ̇(t_n+1)    ≈(I+h/2μ)^-1[(I-h/2μ)θ̇(t_n)+hT(t_n)+h^2/2Ṫ(t_n)],

and 

    θ(t_n+1)≈θ(t_n)+h/2[θ̇(t_n)+θ̇(t_n+1)].

Therefore, θ(t_n+1) and θ̇(t_n+1) can be
computed using θ(t_n), θ̇(t_n), T(t_n),
Ṫ(t_n), but without using unknowns of the t_n+1 step.



 §.§ Concatenating the whole model


With modal transformation and time discretization of the ODEs established,
it now suffices to combine all parts of the piano model together.
The whole computation process consists of two major separate parts:
first do space discretization, then do time discretization. Rooted
in the time-space separation idea of FEM, these two numeric works
do not interfere so computation costs should be acceptable. 



  §.§.§ Space discretization




  
Derive ODEs for each subsystem.

    Mξ̈(t)+Cξ̇(t)+Kξ(t)=f(t), superscripts:a,b,e,f
       Mξ̇(t)+Kξ(t)=f(t), superscripts:c
       θ̈=-μθ̇+T(θ), superscripts:d





  
Solve generalized eigenvalue problems.

Only find the lowest M eigenvalues and corresponding eigenvectors,
and integrate them into Λ, Φ.

    Kϕ_i=λ_iMϕ_i, i=1,...,M, superscripts:a,b,c,e,f





  
Compute modal force transformation matrix (dense).

This is useful for transformation p(t)=Sf(t).

    S=(Φ^HMΦ)^-1Φ^H, superscripts:a,b,c,e,f





  §.§.§ Time discretization




  
Update schemes for modal DOFs

Below lists the previously introduced time discretization shemes to
be later referred to. For each time step n, schemes 1a/2a/3a are
preferred over 1b/2b/3b whenever possible, because they utilize values
of the rhs source term for next time step n+1.


  * Update scheme 1a:

    q(t_n+1)=(I+h/2Λ)^-1[(I-h/2Λ)q(t_n)+h/2p(t_n)+h/2p(t_n+1)]


  * Update scheme 1b:

    q(t_n+1)=(I+h/2Λ)^-1[(I-h/2Λ)q(t_n)+hp(t_n)+h^2/2ṗ(t_n)]


  * Update scheme 2a:

    q̇(t_n+1)     =Z_1^-1[Z_0q̇(t_n)-hΛq(t_n)+h/2p(t_n)+h/2p(t_n+1)]
    q(t_n+1)     =q(t_n)+h/2[q̇(t_n)+q̇(t_n+1)]


  * Update scheme 2b:

    q̇(t_n+1)     =Z_1^-1[Z_0q̇(t_n)-hΛq(t_n)+hp(t_n)+h^2/2ṗ(t_n)]
    q(t_n+1)     =q(t_n)+h/2[q̇(t_n)+q̇(t_n+1)]


  * Update scheme 3a: 

    θ̇(t_n+1)     =(I+h/2μ)^-1[(I-h/2μ)θ̇(t_n)+h/2T(t_n)+h/2T(t_n+1)]
    θ(t_n+1)     =θ(t_n)+h/2[θ̇(t_n)+θ̇(t_n+1)]


  * Update scheme 3b: 

    θ̇(t_n+1)     =(I+h/2μ)^-1[(I-h/2μ)θ̇(t_n)+hT(t_n)+h^2/2Ṫ(t_n)]
    θ(t_n+1)     =θ(t_n)+h/2[θ̇(t_n)+θ̇(t_n+1)]





  
Initialization


  * Discrete time interval h, number of time steps n_1.

  * Hammer shank: initial angle and angular velocity θ^(d')(t_0),
θ̇^(d')(t_0), from which the initial torque T^(d')(t_0)
and its derivative Ṫ^(d')(t_0) (gravity contribution
only) can be computed.

  * Hammer felt & string & soundboard & air & room barraiers: initial
modal DOFs and modal forces 

    q^(*)(t_0)=q̇^(*)(t_0)=p^(*)(t_0)=ṗ^(*)(t_0)=0, *=a,b,c,e”,f





  
At each time step n=0,1,2,...,n_1


  * Compute the next rotation of hammer shank using scheme 3b:

    θ^(d')(t_n),θ̇^(d')(t_n),T^(d')(t_n),Ṫ^(d')(t_n)→θ^(d')(t_n+1),θ̇^(d')(t_n+1)


  * Compute the current modal forces of hammer felt (non-zero boundary
condition) using (<ref>):

    θ^(d')(t_n),θ̇^(d')(t_n),q^(a)(t_n),q̇^(a)(t_n)→p^(e”)(t_n),ṗ^(e”)(t_n)

Note: p^(e”)(t_n) relates to the displacement
of felt at the contact points with string.

  * Compute the next modal DOFs of hammer felt using scheme 2b: 

    q^(e”)(t_n),q̇^(e”)(t_n),p^(e”)(t_n),ṗ^(e”)(t_n)→q^(e”)(t_n+1),q̇^(e”)(t_n+1)


  * Compute the next torque of hammer shank using (<ref>):

    θ^(d')(t_n+1),θ̇^(d')(t_n+1),q^(e”)(t_n+1),q̇^(e”)(t_n+1)→ T^(d')(t_n+1),Ṫ^(d')(t_n+1)


  * Compute the next modal forces of string (non-conservative force) using
(<ref>):

    q^(e”)(t_n+1)→p_1^(a)(t_n+1)

Note: p_1^(a) relates to hammer felt force, p_2^(a)
relates to bridge point displacement.

  * Compute the next modal DOFs of string using scheme 2a for p_1^(a)
and scheme 2b for p_2^(a): 

    q^(a)(t_n),q̇^(a)(t_n),p_1^(a)(t_n),p_1^(a)(t_n+1),p_2^(a)(t_n),ṗ_2^(a)(t_n)→q^(a)(t_n+1),q̇^(a)(t_n+1)


  * Compute the next modal forces of soundboard (non-conservative force)
using (<ref>):

    q^(a)(t_n+1)→p_1^(b)(t_n+1)

Note: p_1^(b) relates to string force at the bridge,
p_2^(b) relates to air pressure on the soundboard.

  * Compute the next modal DOFs of soundboard using scheme 2a for p_1^(b)
and scheme 2b for p_2^(b): 

    q^(b)(t_n),q̇^(b)(t_n),p_1^(b)(t_n),p_1^(b)(t_n+1),p_2^(b)(t_n),ṗ_2^(b)(t_n)→q^(b)(t_n+1),q̇^(b)(t_n+1)


  * Compute the next modal forces of air (non-zero boundary condition)
using (<ref>):

    q^(b)(t_n+1)→p_1^(c)(t_n+1)

Note: p_1^(c) relates to displacements at the
interface with soundboard, p_2^(c) relates to
displacements at interface with room barriers.

  * Compute the next modal DOFs of air using scheme 1a for p_1^(c)
and scheme 1b for p_2^(c): 

    q^(c)(t_n),p_1^(c)(t_n),p_1^(c)(t_n+1),p_2^(c)(t_n),ṗ_2^(c)(t_n)→q^(c)(t_n+1)


  * Compute the next modal forces of soundboard (non-conservative force)
using (<ref>):

    q^(c)(t_n+1)→p_2^(b)(t_n+1),ṗ_2^(b)(t_n+1)


  * Compute the next modal forces of room barriers (non-conservative force)
using (<ref>):

    q^(c)(t_n+1)→p^(f)(t_n+1)

Note: p^(e”)(t_n) relates to air pressure on
the room barriers.

  * Compute the next modal DOFs of room barriers using scheme 2a: 

    q^(f)(t_n),q̇^(f)(t_n),p^(f)(t_n),p^(f)(t_n+1)→q^(f)(t_n+1),q̇^(f)(t_n+1)


  * Compute the next modal forces of air (non-zero boundary condition)
using (<ref>):

    q^(f)(t_n+1),q̇^(f)(t_n+1)→p_2^(c)(t_n+1),ṗ_2^(c)(t_n+1)


  * If higher accuracy is desired, repeat the above steps 1 to 14 several
times. Note that different from the first iteration, the subsequent
iterations always use schemes 1a, 2a, 3a and do not need to use schemes
1b, 2b, 3b. If no more repetition is needed, end the current time
step and move on to the next time step. 




  
Audio output

Compute the acoustic pressure signals at certaining listening positions
using the stored modal DOFs q^(c)(t_n) (n=0,...,n_1)
and modal superposition. These signals are the final output digital
audio of piano simulation model.



§ CONCLUSION


This paper presented a detailed physical model for simulating acoustic
piano sounds. For solid parts of the piano system, viz. strings, soundboard,
room barriers, hammer felt, a 3D prestressed elasticity model is generally
applied. For fluid parts of the piano system, viz. sound radiation
in the air, conservation of mass and Navier-Stokes equation is applied.
For coupling between different subsystems of the piano, mechanisms
of surface force transmission and displacement/velocity continuity
are considered. For numeric simulation, modal superposition and explicit
time discretization schemes are utilized. Despite the complexity of
this whole piano model, we have paid efforts to a straightforward
presentation based more on system ODEs transformed from strong PDEs,
as well as a time domain simulation scheme balancing efficiency and
accuracy.

Below discusses the current study's limitations and our plans or recommendations
for future research.


  * Waiting for numeric simulation results. Due to the complication of
our piano model, we choose to write down theoretic models first as
a guiding framework. Our next step involves implementing the computation
procedures using high-performance, expressive and well-structured
programming languages like Rust, as well as performing result analysis
using convenient and ecologically rich programming languages like
Python. Facing some unknown uncertainties in practice, our model needs
to be further tested and improved.

  * In hammer felt-string coupling and string-soundboard coupling, the
contact was treated as occuring at a point rather than an area. This
simplification would probably lead to a loss in realism. But unlike
for soundboard-air coupling, the interface region is relatively small
for FEM meshes, thus special efforts may be required to realize surface
force transmission and displacement/velocity continuity along small
regions.

  * The string-soundboard coupling mechanism was assumed of fixation but
not collision nature. A coupling model similar to the nonlinear hammer-string
interaction may be more suitable, considering that a string seems
to actually be supported between two distanced bridge pins. Here a
more general framework for collisions in musical instruments may be
applied <cit.> that may better be implemented
through implicit time schemes. It also remains to discover how the
relative positions of the two pins affect the vibrations of the string
and soundboard.

  * Formulation of the piano model lacks some energy perspectives, since
the authors are currently unfamiliar with analyzing energy in dissipative
systems. Besides, a Lagrangian or Hamiltonian formulation of the model
may be a good complement to our Newtonian formulation, unveiling the
conservation and dynamics of kinetic and potential energy under the
least action and virtual work principles. 

  * Some observed phenomena of acoustic pianos are still missing their
representations in our model. For example, strings not striked by
the hammer may also vibrate as long as the sustain pedal is pressed,
which is often called sympathetic resonance. An explanation for this
is that the striked string transmits its vibration to other strings
through the variation of air pressure, which can actually lead to
a string-air coupling model. Also, the transmission of piano player's
key action into hammer shank movement requires further investigation.

  * It remains to discover the relation between physical parameters of
the piano model and the objective (waveforms, spectrums) and subjective
(listener feel) aspects of the final ouput sound, so that these parameters
can be tuned to approach realism or even obtain new sounds without
a real world couterpart. 

Finally, we hope this study could contribute to the understanding
of the vibroacoustics of a piano, and to the innovation of digital
musical instruments with desired characteristics of sound.



plain







§ PRESTRAIN AND PRESTRESS


Consider 3 configurations of a material: natural → initial → current
<cit.>. In the natural configuration, no prestrain
and prestress is present. In the initial configuration, prestrain
and prestress exist. In the current configuration, strain and stress
emerge in addition to prestrain and prestress. Under the same xyz
coordinate system, define ℝ^3 vectors x̅,
x, x̂ as the coordinates of the
same material particle in the natural, initial and current configurations
respectively. The following relation holds: 

    x    =x̅+u̅(x̅),
    x̂    =x+u(x),

where u̅ and u are the displacement
vectors from natural to initial and from initial to current configurations
respectively. Denote deformation gradients V̂=∂x̂/∂x̅,
V=∂x̂/∂x
and V̅=∂x/∂x̅;
denote Jacobian matrices J=∇u, J̅=∇u̅,.
We then have 

    V=I+J, V̅=I+J̅,
       V̂=VV̅=(I+J)(I+J̅).

The Green-Lagrange strain tensor of current configuration with respect
to natural configuration is 

    Ê    =1/2(V^⊤V-I)
        =1/2(J̅+J̅^⊤+J̅^⊤J̅)+1/2(I+J̅^⊤)(J+J^⊤+J^⊤J)(I+J̅)
       ≈E=ϵ+E̅+(J̅^⊤ϵ+ϵJ̅+J̅^⊤ϵJ̅),

where E̅=1/2(J̅+J̅^⊤+J̅^⊤J̅)
is the prestrain tensor and ϵ=1/2(J+J^⊤)
is the strain tensor approximated as the well-known engineering strain;
the second-order 1/2J^⊤J
term is discarded for the sake of linearization. Assuming J̅={J̅_ij}
is symmetric, viz. no rigid body rotation, then we have eigen decomposition
J̅=Udiag(λ_1,λ_2,λ_3)U^⊤
and 

    J̅+1/2J̅^2    =E̅
    U[[ λ_1+1/2λ_1^2;              λ_2+1/2λ_2^2;                           λ_3+1/2λ_3^2 ]]U^⊤    =U[[ λ_1';      λ_2';           λ_3' ]]U^⊤,

from which λ_i=4√(λ_i'+4)-8 is the solution
and J̅ can be computed per eigen decomposition.
The uniqueness of this solution stems from that the principal stretches
λ_i'≥-1 should hold for normal cases, viz. no negative
compression, and λ_i should be close to λ_i'
to be consistent with the case of ignoring 1/2J̅^2.
If the prestrain is not large enough to induce non-negligible geometric
nonlinearity, 1/2J̅^⊤J̅
can be disregarded and J̅=E̅
simply holds. But generally for string instruments like piano, prestrain
may be large and even dominate the post-strain, thus we should cover
the geometric nonlinearity of prestrain. This, however, would not
necessarily make (<ref>) nonlinear with respect
to J, because strain deformation is normally small
enough to make 1/2J^⊤J negligible.
Consequently, with J̅ and E̅
known, the prestrain model (<ref>) is reasonably
linear with respect to the unknowns. Converting it into vector form,
we find 

    Ê   ≈E=ϵ+E̅+Ψϵ,
    Ψ    =[[                 J̅_11(J̅_11+2)                        J̅_12^2                        J̅_13^2                 J̅_12(J̅_11+1)                 J̅_13(J̅_11+1)                     J̅_12J̅_13;                        J̅_12^2                 J̅_22(J̅_22+2)                        J̅_23^2                 J̅_12(J̅_22+1)                     J̅_12J̅_23                 J̅_23(J̅_22+1);                        J̅_13^2                        J̅_23^2                 J̅_33(J̅_33+2)                     J̅_13J̅_23                 J̅_13(J̅_33+1)                 J̅_23(J̅_33+1);                2J̅_12(J̅_11+1)                2J̅_12(J̅_22+1)                    2J̅_13J̅_23 J̅_11J̅_22+J̅_11+J̅_12^2+J̅_22    J̅_11J̅_23+J̅_12J̅_13+J̅_23    J̅_12J̅_23+J̅_13J̅_22+J̅_13;                2J̅_13(J̅_11+1)                    2J̅_12J̅_23                2J̅_13(J̅_33+1)    J̅_11J̅_23+J̅_12J̅_13+J̅_23 J̅_11J̅_33+J̅_11+J̅_13^2+J̅_33    J̅_12J̅_33+J̅_12+J̅_13J̅_23;                    2J̅_12J̅_13                2J̅_23(J̅_22+1)                2J̅_23(J̅_33+1)    J̅_12J̅_23+J̅_13J̅_22+J̅_13    J̅_12J̅_33+J̅_12+J̅_13J̅_23 J̅_22J̅_33+J̅_22+J̅_23^2+J̅_33 ]]

where all the vectors are defined similar to (<ref>).
It is now straight forward that the prestrain functions a linear transformation
(addition and scaling) of the strain at a first approximation. However,
as in (<ref>) ϵ̅_ij(x̅)
is expressed in natural coordinates, converting it into initial coordinate
representations to align with ϵ
would be preferrable. In order for this, define prestress vector σ̅(x̅)
and constitutive matrix D̅ similar to (<ref>),
then E̅=D̅^-1σ̅.
From (<ref>), an inverse mapping x̅=f(x)
should exist. If we let T(x)=σ̅(f(x)),
the vectorized tension field from (<ref>), and
subsititute E̅=D̅^-1T(x)
into (<ref>), then the total strain Ê
can be expressed in only the initial configurations. The advantage
of this is not only avoiding finding the intricate natural coordinates,
but also making it straightfoward to satisfy static equilibrium by
condition ∇·T(x)=0. 

Given the total strain vector, the total stress vector is 

    σ̂(x)=DÊ≈DE=σ+T+τ, τ=DSϵ,

where τ is the dynamic prestress
vector as a secondary contribution of prestress to the total stress;
all vectors here can be converted to corresponding symmetric matrices.
The total stress here is the Cauchy stress, also known as the true
stress, which represents surface forces in the current configuration.
This may cause inconsistencies as σ̂
is still expressed in the initial configuration in (<ref>),
and a conversion to the current configuration may be inconvenient
to apply. Nevertheless, the differences from this should be negligible
under small deformations. It is now clear that the effect prestress
can be seen as bringing a modified non-symmetric constitutive matrix
D(I+S).

In the presence of large deformations and geometric nonlinearity,
it may be beneficial to use the second Piola-Kirchhoff stress tensor
(PK2) (V̂)V̂^-1σ̂V̂^-⊤
<cit.>. This alternative symmetric stress tensor
is derived from the energy conjugate to the Green-Lagrange strain,
capable of representing stress in a pre-deformation configuration,
but seems less physically clear and interpretable compared to the
Cauchy stress tensor. If using PK2, we would seek for the initial
configuration representation as it is used throughout the whole piano
model, thereby the PK2 stress is (V)V^-1σ̂V^-⊤,
which can be first-order or second-order approximated to reduce complexity.
Notice, however, that even when the prestrain is large, the total
stress is already represented in the initial configuration which exists
after deformation from the natural state, eliminating the need to
go back to the natural configuration. And also considering that the
“post” strain would normally be not too large, the difference
between the initial and current configurations should be under an
acceptable level. Therefore, the Cauchy stress tensor is appropriate
for our piano model. As for the choice between engineering strain
(linear) and Green-Lagrange strain (linear + nonlinear) in (<ref>),
it depends on whether prestrain or strain is large enough to induce
non-negligible geometric nonlinearity.


