Popular Models and Numerical Methods for Propagation of Acoustical Shock Waves

Get Complete Project Material File(s) Now! »

Motivation and Objective

Instances involving propagation of acoustical shock waves in complex geometry are nu-merous, sometimes-wanted and sometimes-unwanted. Here are a few examples presented which motivate this thesis project.

Buzz-Saw Noise

As defined by McAlpine et al. [103], it is the noise generated from the turbo-fan engine of an aeroplane when the relative speed of the inlet flow impinging on the fan blades is supersonic. The pressure field associated to a supersonic ducted fan, in a direction normal to the shock fronts, looks like a sawtooth waveform. This is how, this buzzing noise gets the name Buzz-saw noise. It is particularly significant during the take-off and climb, and affects the sound level of the cabin and community. With the increase of commercial air traffic, it becomes important to predict and control these emissions. It has been previously discussed by several authors like Philpot [109], and Hawkings [67].This nonlinear propagation of high amplitude sawtooth wave form inside a duct is classical example of propagation of acoustical shock waves in complex geometry. Here, the complex geometry is due to the inner shape of the turbofan.


Sonic Boom

The introduction of supersonic flights in 1950’s brought into the phenomenon of sonic boom, it is a well-known phenomenon in the field of acoustical shock waves. The pressure disturbance created by the supersonic jet transforms into a N-wave i.e., weak acoustical shock waves. This N-wave is annoying for the population. A detailed discussion on the nature of sonic boom is done in [97, 110, 137]. The interaction of N-wave with topography can lead to diffraction of shock waves and formation of a shadow zone [11, 40]. This interaction brings in the complex geometries as it could be any landscape.

Reflection of Shock Waves

Reflection of acoustical shock waves over a rigid surface are one of the most fundamental phenomenon for shock waves, including acoustical shock waves. The reflection can be broadly classified into two categories namely, regular and irregular (see Ben-Dor [10]). The type of reflection depends on the grazing angle and the strength of the incident shock.
Regular reflection is one which has 2 shock fronts which are the incident and the reflected fronts. It is observed for a sufficiently large grazing angle or a sufficiently weak shock. It is further subdivided into two categories. First case is when reflection obeys the linear Snell-Descartes law of reflection, secondly, when the reflected shock has a curvature and therefore has a varying angle of reflection.
As the criterion for regular reflection is no more satisfied, the point of intersection of the incident and the reflected shock detaches form the surface and gives rise to a third shock. This new shock is called the Mach shock/stem [96]: it connects the merging point of the two shocks and the surface, and is called the triple point. It is important to mention that the slope has a discontinuity at the triple point. Such type of reflection comes under the category of irregular reflection. Premier theory for shock wave reflection was done by von Neumann [133], he called irregular reflection as three-shock theory and regular reflection as two-shock theory.
Colella and Henderson [36] observed numerically and experimentally that for weak shocks there is no triple point, the reflected shock front has a continuous slope along the incident shock and the Mach shock. This happens because the reflected shock breaks down in a band of compressive waves as it approaches the incident shock. They called this new type of reflection as von Neumann reflection. Such weak shock waves exist in acoustics. First numerical observation of nonlinear reflections are done by Sparrow et al.
[121]. Baskar et al. [8] studied the transition in detail theoretically and numerically. They observed the one-shock irregular reflection at almost grazing case where the reflected is not visible as it merges with the incident shock. They call it as weak von Neumann reflec-tion. An experimental validation is done underwater by Marchiano et al. [98]. Karzova et al. [81, 82] studied the interaction of weak shock waves leading to formation of Mach stem in focused beams using optical instruments. Moreover, Pinton et al. [111] simulated the nonlinear reflection of acoustical shear shock waves in soft elastic tissues (involving cubic nonlinearities). This application demonstrates the importance of nonlinear effects near the region of reflection in the propagation of shock waves.


Ultrasound has gained importance for therapeutic applications. Extracorporeal shock wave lithotripsy (ESWL) is used for breaking stones in human body, when they are too big to pass through the urinary tract. It is the most prominent example of therapeutic ultrasound. The first attempts for such a procedure were made in 1950s [88]. It has been successfully implemented since 1980 [24] for fragmenting kidney stones, and later on for gall bladder stones [117].
ESWL involves focusing of high amplitude acoustical shock waves that are generated outside the body and are focused onto a stone within the body. Due to the focusing there is a very high pressure on the stone and significantly lower in the surrounding. The patient is positioned in a way such that the focus of the lithotripter coincides with the stone inside the body, this is achieved through ultrasonic imaging (or other imaging devices). There have been various explanations of the destruction of the stone like compressive failure [23], spalling [48], cavitation [38, 44]. The problems associated with lithotripsy includes Hematuria, renal injury [83, 51], spalling in tissues at the air interfaces such as lung [46] and intestines [64]. Cavitation is also associated to the injury bubble implosion could lead to tissue damage [37].
Dornier HM3 is one of the first and most popular lithotripter in clinical and scientific fraternity [28, 26]. The geometry includes an ellipsoidal geometry with one focus as the source of shock waves and the other focus is made at the stone location inside the body. In other words, half-ellipsoid is outside the body and acts as the mirror to reflect and focus the shock waves at the stone and breaks it (all this is done without any surgery). This clearly demonstrates that it is a well-suited example for propagation of acoustical shock waves in complex geometry.

High Intensity Focused Ultrasound

As mentioned before, use of ultrasound in therapeutic applications is getting importance [5]. High intensity focused ultrasound (HIFU) is used for noninvasive thermal destruction of tumors (see Crum and Hynynen [43], ter Haar [123]), to stop hemorrhage of punctured blood vessels (Vaezy et al. [130]), acoustic characterization Hoff [71], breaking down of microscopic structures Burov et al. [20]. The HIFU devices are constructed using the two-dimensional phased arrays (see Pernot et al. [107], Hand et al. [60]) along a spherical aperture. The ultrasound waves emitted by the transducers are focused on the center of the sphere, which is expected to be a tumor in case of hyperthermic treatments. A detailed discussion on HIFU can be found in [139].
Traditionally, HIFU does not involve shock waves but high intensity continuous waves. Nevertheless, Canney et al. [22] showed that the use of shock waves can improve the heating effects. Indeed, higher frequencies are more readily absorbed and converted to heat than the fundamental frequency. Therefore, the impact of enhanced heating due to acoustical shock waves could be either useful or dangerous and should be properly estimated. It is important to note that this example involves a complex geometry as the medium could have bones and other tissues, and thermal effect could damage them severely. This makes it an interesting case for shock propagation, especially if there are heterogeneities in the domain.
Above examples illustrate different situations involving acoustical shock waves in com-plex geometries. Although, this is not an exhaustive list but it illustrates well the variety of problems which motivates this work and the need of a numerical solver for the propa-gation of acoustical shock waves in complex geometry.

Popular Models and Numerical Methods for Propagation of Acoustical Shock Waves

Numerous works on the propagation of acoustical shock waves using different models have been done since the beginning of numerical computing. In this section a rough survey is done for different models and their associated numerical methods. We choose to present them in an increasing order of complexity starting from the simplest 1D equation to the most general system of equations. Note that, this order corresponds more or less to the historical development of numerical simulation of shock waves.
The simplest equation for propagation of acoustical shock waves is the inviscid Burgers equation [113, 18]. It is a 1D nonlinear advection equation. Starting from a smooth initial waveform, it can take into account the steepening of the waveform until the formation of a discontinuity called the acoustical shock. Once the shock is formed the Burgers equation alone cannot manage the shock as it could lead to multi-valued solution and so the weak shock theory is coupled to provide a physically admissible solution [59, 137]. Beyond 1D problem, it is also used to solve multi-dimensional problems like sonic boom coupled with the technique of ray-tracing. Many numerical methods have been used to solve this equation, details can be found in the textbooks [125, 92, 93, 70, 58]. Note that in this work, we solve the Burgers equation for the development of the method. To assess the quality of the numerical solution, we compare the solution using a Burgers-Hayes quasi-analytical solution developed by Coulouvrat [41] based on so-called Burgers-Hayes method [19, 68].
The next model is the Khokhlova-Zabolotskaya-Kuznetsov (KZK) equation [87] or the KZ equation (KZK without the thermoviscous effects) [141].This is a one-way equation which takes into account the diffraction, nonlinearity and attenuation with a limited angular validity. Indeed, its derivation is based on paraxial approximation of the propa-gation operator. Note that, this equation can be reduced to the Burgers equation if the diffraction is not taken into the account. This model is very useful to simulate propaga-tion of narrow beams in acoustics. The first implementation has been the calculation of the pressure field produced by axisymmetric sources in the near field of a piston completely in frequency domain [1]. This code is the known as the Bergen code. It was later used to investigate the focused beams [61], interaction between finite amplitude beams [124]. Three dimensional codes were developed [80, 21] to investigate the generation of harmonics from a rectangular aperture source. This spectral method was very successful, especially in handling attenuation but was not efficient for strong nonlinearity, where there is generation of higher harmonics (Gibbs phenomenon). There are three other pop-ular approaches based on the fractional step procedure [4]. Since KZK is a one-way wave equation, it models the propagation in the privileged direction. It involves splitting of the physical effects in each spatial advancement step. The same procedure is carried out iteratively. The first of its kind was proposed by Bakhvalov et al. [6]. It solves diffraction and attenuation in frequency domain and nonlinearity in temporal domain. This pseudo-spectral method has also been used to treat the problem of focusing of sonic boom on fold caustics by Marchiano et al. [99, 101] through a generalized KZ equation with the hetero-geneous term proportional to the distance from the caustic. Another method proposed by Lee and Hamilton [90], solves the KZK equation directly in time domain. This code is known as the Texas code. Coulouvrat and co-workers [42, 100] used split-step approach to study the nonlinear Fresnel diffraction and focusing of shock waves. Conclusively, the popularity of this model is due to the fact that its simulation is really fast and efficient but is limited to paraxial approximation.
Several improvements have been proposed to go beyond the parabolic approximation. First of all, the wide-angle approximation [27] is done to extend the angle of validity [56]. Christopher and Parker [25] proposed a method without any angular restriction which relies on the phenomenological way. Recently, Dagrau et al. [45] introduced the HOWARD method which stands for heterogeneous one-way approximation for the reso-lution of diffraction. The numerical resolution is based on the pseudo-spectral approach, diffraction and heterogeneities are solved using spectral methods and nonlinear effects are taken into account using Burgers-Hayes analytical solution [68, 41]. It has been extended to simulate the propagation of shock waves in flows (FLHOWARD [55]). Nevertheless, though these methods have no limitation of angular validity, they are still one-way meth-ods. Consequently, they cannot take into account the effects of back scattering due to heterogeneities or boundaries in complex geometries.
Back scattering effects can be taken into account only using a full-wave approach. The simplest model dealing with propagation in all directions of space and nonlinearity is the Westervelt equation [136, 59]. It consists of a scalar wave equation augmented with a nonlinear term similar to the one in Burgers equation. Note that, in the derivation of this equation the local nonlinear effects are not taken into account (see Chapter 2 for details). Therefore, this is not the most general nonlinear wave equation. The most general nonlinear wave equation in fluid is the Kuznetsov equation [87] which incorpo-rates both local and cumulative nonlinear effects. Nevertheless, the most popular is the Westervelt equation because the local nonlinear effects are expected to be small [1, 79] and from a numerical point of view the remaining nonlinear term is simpler to solve. Different numerical techniques exists: Pinton et al. [112] proposed to solve it using the FDTD, Treeby et al. [128] used k-space method, Verweij et al. [132] used the convolu-tion approach. Note that, all these methods are having limitations in handling complex geometries or steep shocks, although possible approaches are in development (see [127] for an implementation of nonuniform grid in 1D).
As mentioned before, the direct resolution of the Kuznetsov equation is not easily implementable. A first-order system of equations equivalent to Kuznetsov is usually pre-ferred. The pioneering work has been proposed by Sparrow et al. [121] who derived a system based on primitive variables (retaining up to the quadratic terms) and solved it using the FDTD method in Cartesian mesh. They showed the formation of Mach stem using a spherical source over a plane surface. Ginter et al. [57] used similar system in axisymmetric form to investigate nonlinear ultrasound propagation in ideal fluids: it was solved using FDTD approach, in which they are using the DRP (Dispersion relation preserving) scheme. Delpino et al. [47] proposed a very high order finite volume method to simulate the propagation of shock waves induced by explosive source in air. Velasco-Seguar and Rendon [131] recently implemented a low-order finite volume method based on the CLAWPACK codes [92] on graphical processing units, but it requires fine dis-cretization to capture the shock. Few researchers are solving directly the Euler [138, 102] or the Navier-Stokes [3] equations for the propagation of nonlinear waves.
Again all these examples deal with regular geometries. Nevertheless, there are clear advantages of solving the system of first order equations. It is closer to the physics than wave equations (conservation properties). It gives access to all the velocity components, density variations and the pressure. This enables a more detailed study of different phe-nomenon of reflection, refraction, diffraction, attenuation, dispersion, nonlinearity. For instance, the effect of Lagrangian density [1], which is a local effect, can be studied accu-rately. Nevertheless, as it has been outlined, it is difficult to handle complex geometries. A solution is to use a method build on unstructured mesh. To our knowledge, such a method has not yet been developed for the system of nonlinear equations.


Numerical Methods for Complex Geometry and Acoustical Shock Waves

Choice of the Method

In this section, we discuss about the main numerical methods and their ability to propa-gate acoustical shock waves in complex geometries. Generally in nonlinear acoustics there are long propagation distance involved (about 100 wavelengths), for which there is a need of a method with low dispersion and low dissipation. Such attributes are contained in a high-order methods. Therefore, it implies three features: high order schemes for long propagation, handling complex-domains, capturing of nonlinear effects including shock formation, propagation and merging.
The finite difference methods (FDMs or FDTD in acoustics) [35] are the most popular methods for solving the nonlinear partial differential equations as seen in the previous section. Indeed, they are easily implementable. It is relatively easy to get high order dis-cretization in space, which gives the freedom to choose an efficient time-stepping method.
These features make it applicable to variety of problems in nonlinear acoustics. However, despite techniques like curvilinear coordinates or immersed boundary condition, the finite difference methods are ill-equipped for handling complex geometries [52, 129].
The family of FVMs are good in handling the problem of complex geometry as, in these methods, the space is discretized in volumes or cells. In each cell, the numerical computations are purely local and the fluxes are computed with the neighboring cells. Higher-order spatial accuracy in finite volume methods involves re-construction of cell averages. This creates an expanded numerical stencil, which drastically impacts the itera-tive algorithm and also complicates implementation of boundary condition. Nevertheless, FVMs are the most popular methods for hyperbolic problems, with second-order accurate methods being most frequently used [93, 125].
On the other hand in the FEMs [73, 142, 122], a spectral solution is constructed using a globally defined basis and with the same test functions. This gives a implicit semi-discrete form and the mass matrix is required to be inverted. Here the problem is the large global mass matrix which requires large memory. Moreover, it could also lead to instabilities [69]. Such methods are the best choice for problems like heat equation but not for wave propagation problems.
The Discontinuous Galerkin Method (DGM) is a kind of hybrid of the FEM and the FVM. It is capable of handling complex geometries thanks to unstructured mesh. DGM preserves the spectral nature of the solution within one element as in FEMs based on basis and test functions, and can have high order representation. But, it satisfies the equations locally within each element this attribute resembles the FVMs. This gives DGM the ability of local (within a element) high-order accuracy, wherever needed. Therefore, it happens to be an appealing choice.
The DGM was first proposed by Reed and Hill [115] for solving a steady-state neutron transport equation, with its analysis given by Lesaint and Raviart [91]. At present the DGM is widely applied to many areas [69]. In acoustics, it has been used mainly for linear acoustics [85], aeroacoustics [126, 53, 54], propagation at the interface between moving media and isotropic solids Luca et al. [94, 95], and nonlinear acoustics in solids [15]. To our knowledge, DGM has not been used for propagation of acoustical shock waves. Indeed, when an acoustical shock appears the method does not capture it properly by itself.

Table of contents :

1 Introduction
1.1 Motivation and Objective
1.2 Popular Models and Numerical Methods for Propagation of Acoustical Shock Waves
1.3 Numerical Methods for Complex Geometry and Acoustical Shock Waves
1.3.1 Choice of the Method
1.3.2 Shock Management
1.4 Outline of the Manuscript
2 Equations of Propagation in Nonlinear Acoustics
2.1 Conservation Laws
2.2 Equations for Nonlinear Acoustics
2.3 Dimensionless Formulation of the System of Equations
2.3.1 Characteristic Parameters and Variables
2.4 Summary
2.5 Comparison with other Equations of Nonlinear Acoustics
2.5.1 Conservative to Primitive form
2.5.2 Kuznetsov Equation
2.5.3 Westervelt Equation
2.5.4 KZ Equation
2.5.5 Inviscid Burgers Equation
2.6 Conclusions
3 Discontinuous Galerkin Method
3.1 Nodal Discontinuous Galerkin Method in 1D
3.1.1 Weak Formulation
3.1.2 Computations in Reference Element
3.1.3 Assembling
3.2 Nodal Discontinuous Galerkin Method in 2D
3.2.1 Weak Formulation
3.2.2 Computations in Reference Element
3.2.3 Assembling
3.3 Brief Review on GPU Implementation
3.4 Time Discretization
3.5 Application 1D: Advection Equation
3.6 Conclusions
4 Shock Management in One-Dimension
4.1 Illustration
4.2 Slope Limiters
4.2.1 Slope Limiter: Cockburn
4.2.2 Slope Limiter: Biswas
4.2.3 Slope Limiter: Burbeau
4.2.4 Numerical Experiment
4.3 Method of Global Artificial Viscosity
4.3.1 Local Discontinuous Galerkin Method in 1D
4.3.2 Numerical Experiment
4.4 Element Centered Smooth Artificial Viscosity
4.4.1 Shock Sensor
4.4.2 Smooth Artificial Viscosity
4.4.3 Implementation Issues
4.5 Validation
4.5.1 Inverted Sine-period to N-wave
4.5.2 Sine-period to Sawtooth
4.5.3 N-wave
4.5.4 Sawtooth
4.5.5 Multiple Shocks
4.6 Conclusions
5 Shock Management in Two-Dimensions
5.1 Equations of Nonlinear Acoustics
5.2 Convective-Diffusive System for Nonlinear Acoustics
5.3 Local Discontinuous Galerkin Implementation
5.3.1 Weak Formulation
5.3.2 Numerical Fluxes
5.3.3 Nodal Approximation
5.3.4 Assembling
5.4 Element Centered Smooth Artificial Viscosity
5.4.1 Shock Sensor
5.4.2 Smooth Artificial Viscosity
5.5 Numerical Explanation of the Shock Sensor
5.5.1 First-Order Contribution to the Shock Sensor
5.5.2 Highest-Order Contribution to the Shock Sensor
5.6 Implementation Issues and Validation
5.7 Conclusions
6 Applications
6.1 Reflection of Acoustical Shock Waves
6.1.1 Numerical Experiments
6.1.2 Results and Discussion
6.2 Focusing of continuous (shock) waves: application to HIFU
6.2.1 Mesh Renement Based on ECSAV
6.2.2 Low resolution simulation
6.2.3 Local high resolution mesh
6.2.4 Focusing in a homogeneous medium
6.2.5 Intensity near the focus
6.2.6 Focusing in a medium with an obstacle
6.2.7 Conclusions
7 Conclusions and Perspectives


Related Posts