In this report, we present energy quadratization techniques for a
Hamiltonian system of nonlinear wave equations formulated at order 2 in time.
On a generic system, we present the so-called Invariant Energy Quadratization (IEQ)
and Scalar Auxiliary Variable (SAV) methods as well as their energy conservation
properties and discretization strategies in space and time. Unlike the iterative
techniques commonly used for nonlinear systems to guarantee certain invariances,
these two methods lead to algorithms whose complexity is known in advance and rely
on the simple inversion of a linear system at each time step. In spite of an
unconditional stability and an attractive complexity, the literature mentions
problematic application cases with an uncontrolled accuracy. The numerical
properties (stability, consistency and uniform convergence in time with respect
to the CFL) of schemes obtained by hybridization between theta scheme and
quadratization are studied for two classes of nonlinear terms: a nonlinearity
concerning the solution field, and a nonlinearity concerning its gradient. These
results are then applied to a geometrically exact nonlinear piano string for
which numerical results are presented. The influence of the discretization
parameters is studied and related to the theoretical results. The choices for
the best accuracy or computational cost are illustrated. Some parameters can
induce the space-time non-convergence of the schemes for a nonlinearity on the
gradient, as it is the case for the piano string.
preprint
Stability and space/time convergence of Störmer-Verlet time integration of the mixed formulation of linear wave equations
This work focuses on the mixed formulation of linear wave equations. It provides a proof of stability and convergence of time discretisation of a semi discrete linear wave equation in mixed form with Störmer-Verlet time integration, that is uniform as the time step reaches its largest allowed value for stability (Courant-Friedrich-Levy condition), contrary to the proofs recalled here from the literature.
RR
Space time convergence of implicit discretization strategies for the mixed formulation of linear wave equations
This report presents two parametric implicit strategies for discretizing the mixed formulation of linear wave equations, based for the first on the Störmer-Verlet scheme and for the second on the Crank-Nicolson scheme. If the initial conditions are well chosen, these two strategies are equivalent, and equivalence formulas allow the variables of one scheme to be reconstructed from the other. These two schemes are also equivalent to the đ
Ξ
-scheme discretization of the wave equation formulated at order 2. However, the first scheme has a lower complexity than the second one. Numerical results illustrate these equivalences and these efficiency comparisons.
2022
JSV
Transmission line coefficients for viscothermal acoustics in conical tubes
Viscothermal acoustic propagation in gases contained in rigid straight or conical tubes is considered. Under the assumption that the wavelength is much larger than both the boundary layer thickness and the tube radius, pressure and flow are shown to be solutions of a pair of coupled 1D differential equations, formulated as transmission line equations involving complex loss coefficients. The derivation of these loss coefficients, which is usually accomplished in cylinders, is generalized here to conical geometries. In the well-kown case of circular cylinders, the ZwikkerâKosten (ZK) theory is recovered. For circular cones, the expression of the loss coefficients is derived. It involves complex-order spherical harmonics, instead of Bessel functions for circular cylinders, and makes the hydraulic radius appear as a natural relevant geometrical parameter. We show that replacing the classical radius by the hydraulic radius in the ZK theory provides an affordable and accurate approximation of the analytic model derived for cones. The proposed formulas are used to compute the input impedance of a cone, and compared with a 3D reference. In an ideal setting, using the spherical harmonics or the hydraulic radius in the 1D method accurately approximates the full 3D method, and allows to increase accuracy by approximately two orders of magnitude compared to the ZK theory.
JSV
Direct computation of modal parameters for musical wind instruments
Parameters of modal expressions for the input impedance of wind instruments can be computed from the Telegrapher equations with radiating boundary conditions at the bell. One dimensional finite elements are used for spatial discretization. Provided that the models can be put under a specified form, the modal parameters can be computed directly by solving a generalized eigenvalue problem. Viscothermal effects can also be accounted for, as well as open or closed toneholes. Modal shapes can be visualized along the instrument. Parameters of modal expressions for the input admittance of flute-like instruments can also be computed.
2021
JASA
Dissipative time-domain one-dimensional model for viscothermal acoustic propagation in wind instruments.
Accurate modeling of the acoustic propagation in tubes of varying cross section in musical acoustics must include the effects of the viscous and thermal boundary layers. Models of viscothermal losses are classically written in the frequency domain. An approximate time-domain model is proposed in which all of the physical parameters of the instrument, the bore shape or the wave celerity, are explicit coefficients. The model depends on absolute tabulated constants, which only reflect that the pipe is axisymmetric. It can be understood as a telegrapherâs equations augmented by an adjustable number of auxiliary unknowns. A global energy is dissipated. A time discretization based on variational approximation is proposed along with numerical experiments and comparisons with other models.
@article{thibault:hal-03328715,title={{Dissipative time-domain one-dimensional model for viscothermal acoustic propagation in wind instruments.}},author={Thibault, Alexis and Chabassier, Juliette},url={https://hal.archives-ouvertes.fr/hal-03328715},note={Copyright 2021 Acoustical Society of America. This article may be downloaded for personal use only. Any other use requires prior permission of the author and the Acoustical Society of America.The following article appeared in Thibault, A., \\& Chabassier, J. (2021). Dissipative time-domain one-dimensional model for viscothermal acoustic propagation in wind instruments. The Journal of the Acoustical Society of America, 150(2), 1165-1175 and may be found at https://doi.org/10.1121/10.0005537 .},journal={{Journal of the Acoustical Society of America}},publisher={{Acoustical Society of America}},volume={150},number={2},pages={1165-1175},year={2021},month=aug,doi={10.1121/10.0005537},pdf={https://hal.archives-ouvertes.fr/hal-03328715/file/Manuscript%20%281%29.pdf},hal_id={hal-03328715},hal_version={v1},bibtex_show={true},html={https://doi.org/10.1121/10.0005537},abbr={JASA}}
M2AN
Construction and convergence analysis of conservative second order local time discretisation for linear wave equations
In this work we present and analyse a time discretisation strategy for linear wave equations t hat aims at using locally in space the most adapted time discretisation among a family of implicit or explicit centered second order schemes. The proposed family of schemes is adapted to domain decomposition methods such as the mortar element method. They correspond in that case to local implicit schemes and to local time stepping. We show that, if some regularity properties of the solution are satisfied and if the time step verifies a stability condition, then the family of proposed time discretisations provides, in a strong norm, second order space-time convergence. Finally, we provide 1D and 2D numerical illustrations that confirm the obtained theoretical results and we compare our approach on 1D test cases to other existing local time stepping strategies for wave equations.
AA
Full waveform inversion for bore reconstruction of woodwind-like instruments
The internal geometry of a wind instrument can be estimated from acoustic measurements. For woodwind instruments, this involves characterizing the inner shape (bore) but also the side holes (dimensions and location). In this study, the geometric parameters are recovered by a gradient-based optimization process, which minimizes the deviation between simulated and measured linear acoustic responses of the resonator for several fingerings through an observable function. The acoustic fields are computed by solving a linear system resulting from the 1D spectral finite elements spatial discretization of the wave propagation equations (including thermo-viscous effects, radiation and side holes). The âfull waveform inversionâ (FWI) technique exploits the fact that the gradient of the cost function can be computed by solving the same linear system as that of the direct problem but with a different source term. The gradient is computed with better accuracy and less additional cost than with finite-difference. The dependence of the cost function on the choice of the observed quantity, the frequency range and the fingerings used, is first analyzed. Then, the FWI is used to reconstruct, from measured impedances, an elementary instrument with 14 design variables. The results, obtained in about 1 minute on a laptop, are in excellent agreement with the direct geometric measurements.
This paper deals with the construction of a family of fourth order, energy consistent, explicit time discretizations for dissipative linear wave equations. The schemes are obtained by replacing the inversion of a matrix, that comes naturally after using the technique of the Modified Equation on the second order Leap Frog scheme applied to dissipative linear wave equations, by explicit approximations of its inverse. The stability of the schemes are studied using an energy analysis and a convergence analysis is carried out. Numerical results in 1D illustrate the space/time convergence properties of the schemes and their efficiency is compared to more classical time discretizations.
2019
SPM
Model-based digital pianos: from physics to sound synthesis
As a result of their complexity and versatility, pianos are arguably one of the most important instruments in Western music. The size, weight, and price of grand pianos as well as their relatively simple control surface (i.e., the keyboard), have led to the development of digital counterparts that mimic the sound of acoustic pianos as closely as possible. While most commercial digital pianos are based on sample playback, it is also possible to reproduce a pianoâs sound by modeling the physics of the instrument. The process of physical modeling starts with first understanding the physical principles, followed by creating accurate numerical models, and finally, finding numerically optimized signal processing models that allow sound synthesis in real time by excluding inaudible phenomena and adding some perceptually important features by using signal processing tricks. Accurate numerical models can be used by physicists and engineers to understand how the instrument functions or to help piano makers with instrument development. On the other hand, efficient real-time models are geared toward composers and musicians who perform at home or onstage. This article provides an overview of a physics-based piano synthesis beginning with computationally heavy, physically accurate approach followed by a discussion of the approaches that are designed to produce the best possible sound quality for real-time synthesis.
RR
About the transfert matrix method in the context of acoustical wave propagation in wind instruments
This work presents a computation tool for the calculation of wind instrument input impedance in the context of linear planar wave propagation with visco-thermal losses. The originality of the approach lies in the usage of a specific and simple 1D finite element method (FEM). The popular Transfer Matrix Method (TMM) is also recalled and a seamless formulation is proposed which unifies the cases cylinders vs. cones. Visco-thermal losses, which are natural dissipation in the system, are not exactly taken into account by this method when arbitrary shapes are considered. The introduction of an equivalent radius leads to an approximation that we quantify using the FEM method. The equation actually solved by the TMM in this case is exhibited. The accuracy of the two methods (FEM and TMM) and the associated computation times are assessed and compared. Although the TMM is more efficient in lossless cases and for lossy cylinders, the FEM is shown to be more efficient when targeting a specific precision in the realistic case of a lossy trumpet. Some additional features also exhibit the robustness and flexibility of the FEM over the TMM. All the results of this article are computed using the open-source python toolbox OpenWind.
AMC
Equivalent boundary conditions for acoustic media with exponential densities. Application to the atmosphere in helioseismology
We present equivalent boundary conditions and asymptotic models for the solution
of a transmission problem set in a domain which represents the sun and its atmosphere.
This problem models the propagation of an acoustic wave in time-harmonic regime.
The specific non-standard feature of this problem lies in the presence of a small
parameter ÎŽwhich represents the inverse rate of the exponential decay of the density
in the atmosphere. This problem is well suited for the notion of equivalent conditions
and the effect of the atmosphere on the sun is as a first approximation local.
This approach leads to solve only equations set in the sun. We derive rigorously
equivalent conditions up to the fourth order of approximation with respect to ÎŽfor the exact solution. The construction of equivalent conditions is based on a multiscale
expansion in power series of ÎŽfor u. Numerical simulations illustrate the theoretical results.
Finally we measure the boundary layer phenomenon by introducing a characteristic
length that turns out to depend on the mean curvature of the interface between the subdomains.
2018
RR
Solving time-harmonic Galbrunâs equation with an arbitrary flow. Application to Helioseismology
This work offers some contributions to the numerical study of acoustic waves propagating in the Sun and its atmosphere. The main goal is to provide boundary conditions for outgoing waves in the solar atmosphere where it is assumed that the sound speed is constant and the density decays exponentially with radius. Outgoing waves are governed by a Dirichlet-to-Neumann map which is obtained from the factorization of the Helmholtz equation expressed in spherical coordinates. For the purpose of extending the outgoing wave equation to axisymmetric or 3D cases, different approximations are implemented by using the frequency and/or the angle of incidence as parameters of interest. This results in boundary conditions called atmospheric radiation boundary conditions (ARBC) which are tested in ideal and realistic configurations. These ARBCs deliver accurate results and reduce the computational burden by a factor of two in helioseismology applications.
2017
CRAS
Space/Time convergence analysis of a class of conservative schemes for linear wave equations
This paper concerns the space/time convergence analysis of conservative two-step time discretizations for linear wave equations. Explicit and implicit, second- and fourth-order schemes are considered, while the space discretization is given and satisfies minimal hypotheses. Convergence analysis is done using energy techniques and holds if the time step is upper-bounded by a quantity depending on space discretization parameters. In addition to showing the convergence for recently introduced fourth-order schemes, the novelty of this work consists in the independency of the convergence estimates with respect to the difference between the time step and its greatest admissible value.
AA
Atmospheric radiation boundary conditions for high frequency waves in time-distance helioseismology
The temporal covariance between seismic waves measured at two locations on the solar surface is the fundamental observable in time-distance helioseismology. Above the acoustic cut-off frequency (Â 5.3 mHz), waves are not trapped in the solar interior and the covariance function can be used to probe the upper atmosphere. We wish to implement appropriate radiative boundary conditions for computing the propagation of high-frequency waves in the solar atmosphere. We consider recently developed and published radiative boundary conditions for atmospheres in which sound-speed is constant and density decreases exponentially with radius. We compute the cross-covariance function using a finite element method in spherical geometry and in the frequency domain. The ratio between first- and second-skip amplitudes in the time-distance diagram is used as a diagnostic to compare boundary conditions and to compare with observations. We find that a boundary condition applied 500 km above the photosphere and derived under the approximation of small angles of incidence accurately reproduces the âinfinite atmosphereâ solution for high-frequency waves. When the radiative boundary condition is applied 2 Mm above the photosphere, we find that the choice of atmospheric model affects the time-distance diagram. In particular, the time-distance diagram exhibits double-ridge structure when using a Vernazza Avrett Loeser atmospheric model.
AA
Computational helioseismology in the frequency domain: acoustic waves in axisymmetric solar models with flows
Context. Local helioseismology has so far relied on semi-analytical methods to compute the spatial sensitivity of wave travel times to perturbations in the solar interior. These methods are cumbersome and lack flexibility.
Aims. Here we propose a convenient framework for numerically solving the forward problem of time-distance helioseismology in the frequency domain. The fundamental quantity to be computed is the cross-covariance of the seismic wavefield.
Methods. We choose sources of wave excitation that enable us to relate the cross-covariance of the oscillations to the Greenâs function in a straightforward manner. We illustrate the method by considering the 3D acoustic wave equation in an axisymmetric reference solar model, ignoring the effects of gravity on the waves. The symmetry of the background model around the rotation axis implies that the Greenâs function can be written as a sum of longitudinal Fourier modes, leading to a set of independent 2D problems. We use a high-order finite-element method to solve the 2D wave equation in frequency space. The computation is embarrassingly parallel, with each frequency and each azimuthal order solved independently on a computer cluster.
Results. We compute travel-time sensitivity kernels in spherical geometry for flows, sound speed, and density perturbations under the first Born approximation. Convergence tests show that travel times can be computed with a numerical precision better than one millisecond, as required by the most precise travel-time measurements.
Conclusions. The method presented here is computationally efficient and will be used to interpret travel-time measurements in order to infer, e.g., the large-scale meridional flow in the solar convection zone. It allows the implementation of (full-waveform) iterative inversions, whereby the axisymmetric background model is updated at each iteration.
WaMot
Numerical robustness of single-layer method with Fourier basis for multiple obstacle acoustic scattering in homogeneous media
We investigate efficient methods to simulate the multiple scattering of obstacles in homogeneous media. With a large number of small obstacles on a large domain, optimized pieces of software based on spatial discretization such as Finite Element Method (FEM) or Finite Difference lose their robustness. As an alternative, we work with an integral equation method, which uses single-layer potentials and truncation of Fourier series to describe the approximate scattered field. In the theoretical part of the paper, we describe in detail the linear systems generated by the method for impenetrable obstacles, accompanied by a well-posedness study. For the numerical performance study, we limit ourselves to the case of circular obstacles. We first compare and validate our codes with the highly optimized FEM-based software Montjoie. Secondly, we investigate the efficiency of different solver types (direct and iterative of type GMRES) in solving the dense linear system generated by the method. We observe the robustness of direct solvers over iterative ones for closely-spaced obstacles, and that of GMRES with LowerâUpper Symmetric GaussâSeidel and Symmetric GaussâSeidel preconditioners for far-apart obstacles.
2016
M2AN
Time Domain Simulation of a Piano. Part 2 : Numerical Aspects
This article is the second of a series of two papers devoted to the numerical simulation of the piano. It concerns the numerical aspects of the work, the implementation of a piano code and the presentation of corresponding simulations. The main difficulty is time discretization and stability is achieved via energy methods. Numerical illustrations are provided for a realistic piano and compared to experimental recordings.
RR
High Order Finite Element Method for solving Convected Helmholtz equation in radial and axisymmetric domains. Application to Helioseismology
A family of fourth-order coupled implicitâexplicit time schemes is presented as a special case of fourth-order coupled implicit schemes for linear wave equations. The domain of interest is decomposed into several regions where different fourth-order time discretizations are used, chosen among a family of implicit or explicit fourth-order schemes. The coupling is based on a Lagrangian formulation on the boundaries between the several non-conforming meshes of the regions. A global discrete energy is shown to be preserved and leads to global fourth-order consistency in time. Numerical results in 1D and 2D for the acoustic and elastodynamics equations illustrate the good behavior of the schemes and their potential for the simulation of realistic highly heterogeneous media or strongly refined geometries, for which using everywhere an explicit scheme can be extremely penalizing. Accuracy up to fourth order reduces the numerical dispersion inherent to implicit methods used with a large time step and makes this family of schemes attractive compared with second-order accurate methods.
2014
JSV
Energy based simulation of a Timoshenko beam in non-forced rotation. Application to the flexible piano hammer shank.
A nonlinear model for a vibrating Timoshenko beam in non-forced unknown rotation is derived from the virtual work principle applied to a system of beam with mass at the end. The system represents a piano hammer shank coupled to a hammer head. An energy-based numerical scheme is then provided, obtained by non-classical approaches. A major difficulty for time discretization comes from the nonlinear behavior of the kinetic energy of the system. This new numerical scheme is then coupled to a global energy-preserving numerical solution for the whole piano. The obtained numerical simulations show that the pianistic touch clearly influences the spectrum of the piano sound of equally loud isolated notes. These differences do not come from a possible shock excitation on the structure, or from a changing impact point, or a âlongitudinal rubbing motionâ on the string, since neither of these features is modeled in our study.
The purpose of this study is the time domain modeling of a piano. We aim at explaining the vibratory and acoustical behavior of the piano, by taking into account the main elements that contribute to sound production. The soundboard is modeled as a bidimensional thick, orthotropic, heterogeneous, frequency dependent damped plate, using Reissner Mindlin equations. The vibroacoustics equations allow the soundboard to radiate into the surrounding air, in which we wish to compute the complete acoustical field around the perfectly rigid rim. The soundboard is also coupled to the strings at the bridge, where they form a slight angle from the horizontal plane. Each string is modeled by a one dimensional damped system of equations, taking into account not only the transversal waves excited by the hammer, but also the stiffness thanks to shear waves, as well as the longitudinal waves arising from geometric nonlinearities. The hammer is given an initial velocity that projects it towards a choir of strings, before being repelled. The interacting force is a nonlinear function of the hammer compression. The final piano model is a coupled system of partial differential equations, each of them exhibiting specific difficulties (nonlinear nature of the string system of equations, frequency dependent damping of the soundboard, great number of unknowns required for the acoustic propagation), in addition to couplingsâ inherent difficulties.
A time-domain global modeling of a grand piano is presented. The string model includes internal losses, stiffness, and geometrical nonlinearity. The hammer-string interaction is governed by a nonlinear dissipative compression force. The soundboard is modeled as a dissipative bidimensional orthotropic ReissnerâMindlin plate where the presence of ribs and bridges is treated as local heterogeneities. The coupling between strings and soundboard at the bridge allows the transmission of both transverse and longitudinal waves to the soundboard. The soundboard is coupled to the acoustic field, whereas all other parts of the structure are supposed to be perfectly rigid. The acoustic field is bounded artificially using perfectly matched layers. The discrete form of the equations is based on original energy preserving schemes. Artificial decoupling is achieved, through the use of Schur complements and Lagrange multipliers, so that each variable of the problem can be updated separately at each time step. The capability of the model is highlighted by series of simulations in the low, medium, and high register, and through comparisons with waveforms recorded on a Steinway D piano. Its ability to account for phantom partials and precursors, consecutive to string nonlinearity and inharmonicity, is particularly emphasized.
JCAM
Introduction and study of fourth order theta schemes for linear wave equations
A new class of high order, implicit, three time step schemes for semi-discretized wave equations is introduced and studied. These schemes are constructed using the modified equation approach, generalizing the Ξ-scheme. Their stability properties are investigated via an energy analysis, which enables us to design super-convergent schemes and also optimal stable schemes in terms of consistency errors. Specific numerical algorithms for the fully discrete problem are tested and discussed, showing the efficiency of our approach compared to second order Ξ-schemes.
Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string.
We study the implicit time discretization of piano strings governing equations within the Timoshenko prestressed beam model. Such model features two different waves, namely the flexural and shear waves, that propagate with very different velocities. We present a novel implicit time discretization that reduces the numerical dispersion while allowing the use of a large time step in the numerical computations. After analyzing the continuous system and the two branches of eigenfrequencies associated with the propagating modes, the classical
-scheme is studied. We present complete new proofs of stability using energy-based approaches that provide uniform results with respect to the featured time step. A dispersion analysis confirms that
reduces the numerical dispersion, but yields a severely constrained stability condition for the application considered. Therefore we propose a new
-like scheme, which allows to reduce the numerical dispersion while relaxing this stability condition. Stability proofs are also provided for this new scheme. Theoretical results are illustrated with numerical experiments corresponding to the simulation of a realistic piano string.
This paper considers a general class of nonlinear systems, ânonlinear Hamiltonian systems of wave equationsâ. The first part of our work focuses on the mathematical study of these systems, showing central properties (energy preservation, stability, hyperbolicity, finite propagation velocity, etc.). Space discretization is made in a classical way (variational formulation) and time discretization aims at numerical stability using an energy technique. A definition of âpreserving schemesâ is introduced, and we show that explicit schemes or partially implicit schemes which are preserving according to this definition cannot be built unless the model is trivial. A general energy preserving second order accurate fully implicit scheme is built for any continuous system that fits the nonlinear Hamiltonian systems of wave equations class. The problem of the vibration of a piano string is taken as an example. Nonlinear coupling between longitudinal and transversal modes is modeled in the âgeometrically exact modelâ, or approximations of this model. Numerical results are presented.