Nonlinear Acoustics — Modeling of the 1D Westervelt Equation

Although linear acoustics does an outstanding job of explaining most acoustical phenomena, there are an increasing number of applications that require abandoning the small-signal assumption and applying the wave equation for finite-amplitude sound or nonlinear acoustics. For example, the effects of finite amplitude propagation can be seen in almost every medical use of ultrasound. The use of high intensity ultrasound to induce tissue heating for cancer therapy or bleeding control, and the use of shock waves in extracorporeal and endoscopic lithotripsy of kidney stones and gallstones provide additional examples. Even in the field of diagnostic ultrasound, nonlinear acoustics is of interest because the higher harmonics emerged during wave propagation can be exploited to produce a better image quality.

The nonlinear phenomenon of wave distortion is both a local effect and a cumulative effect due to variation of propagation speed over the waveform, which causes distortion that accumulates with propagation distance. The local effect is usually small compared to cumulative distortion and can be neglected once the propagation distance becomes much greater than a wavelength; see Ref. 1. Therefore, a transient analysis is necessary to model the cumulative distortion along with the wave propagation.

This example shows how to model nonlinear propagation of finite-amplitude acoustic waves in fluids using the Pressure Acoustics, Transient physics interface of the Acoustics Module. The nonlinear effects are taken into account by adding the Nonlinear Acoustics (Westervelt) domain feature. Thus the linear wave equation transforms into the Westervelt equation which is an approximation of the full 2nd-order wave equation when cumulative nonlinear effects dominate local nonlinear effects; see Ref. 1. The current model simulates a finite-amplitude wave propagating in 1D along an interval greater than the shock formation distance. The computed numerical solution is compared to an analytical solution before and after shock formation.

Model Definition

The full Westervelt equation reads

(1)

where p is the total acoustic pressure, ρ0 and c0 are the density and the speed of sound, β = 1 + B/2A is the coefficient of nonlinearity, and δ is the sound diffusivity (see Ref. 2). This is the equation solved when the Nonlinear Acoustics (Westervelt) domain feature is added and the is selected as is selected on the main Transient Pressure Acoustics Model node.

For a piston vibrating in a 1D tube with the velocity u(t) = u0sinωt, the so-called shock formation distance is

There are two classical analytical solutions to Equation 1 available in this case (see Ref. 1). Both of them have the form of a series

(2)

where σ = x/xsh is the dimensionless spatial coordinate and p0 = ρ0c0u0 is the source pressure amplitude.

The first one is known as the Fubini solution and is valid for the distances up to xsh (σ ≤ 1). The harmonic amplitudes Bn are defined as

with Jn being the Bessel function of the first kind of order n. The second one is known as the Fay solution and it is applied for σ ≥ 3.5. The harmonic amplitudes are

with Γ = 2βu0/ωδ is the Goldberg number, which is a measure of the strength of nonlinearity relative to that of dissipation (see Ref. 1). A solution provided by Blackstock exists in the transition region for 1 ≤ σ ≤ 3.5, but this solution is not considered here.

Table 1 shows the material properties of water and some critical parameters used in the simulation. Given those numbers, the shock formation distance for the current model is about 0.1 m. The simulation domain is the interval 0 ≤ x ≤ 4.5xsh, as shown in Figure 1. A sinusoidal pressure source of amplitude p0 is applied at x = 0, and the Plane Wave Radiation boundary condition is applied at x = 4.5xsh to model the propagation of the wave without reflections. The diffusivity of sound present in the model is due to viscosity and is defined as δ = 4μ/3ρ0. It corresponds to the acoustic attenuation α = 8.1·10-5 Np/m.

Table 1: Material properties and some critical model parameters.
Name	Value	Description
ρ0	999.6 kg/m3	Density at 20oC and 1 atm
c0	1481.44 m/s	Speed of sound at 20oC and 1 atm
μ0	1.0016·10-3 Pa·s	Viscosity at 20oC and 1 atm
β	10	Coefficient of nonlinearity
f0	0.1 MHz	Source driving frequency
p0	5 MPa	Source pressure amplitude
N0	8	Number of harmonics to resolve

Figure 1: 1D geometry of the model with the source located at x = 0.

Results and Discussion

The increasing distortion of the waveform when the wave travels away from the acoustic source is illustrated in Figure 2. The nonlinear effects are apparent by comparing the solution with the linear analytical solution. Initially, the waveform distortion is caused by the dependence of propagation speed on the pressure or particle velocity. The peaks of the wave travel faster than the troughs, and the waveform gradually turns into a sawtooth-like shape clearly seen near the right end of the interval. After the shock formation (the red vertical line in Figure 2) the sound dissipation grows. This is due to the intrinsic frequency dependency of the attenuation. It is proportional to the frequency squared.

The distortion of the waveform results in the generation of higher harmonic components. The further the wave travels, the more energy is transferred to the higher harmonic components from the fundamental frequency of the harmonic source signal. This effect is demonstrated in the plots of Figure 3 through Figure 5. The plots compare the model solution (blue lines) with the analytical solutions (green line) for both waveform and frequency spectrum at x = 0.5xsh, xsh, and 3.5xsh. The Fourier transforms used to determine the spectrum are applied for the time interval Δt = 5T0 after the wave arrives at those locations.

At x = 0.5xsh in Figure 3, the 2nd and the 3rd harmonic components start to show up; at x = xsh in Figure 4, more than ten harmonics appear in the frequency spectrum. Their contribution grows with the distance as seen in Figure 5 at x = 3.5xsh. These results clearly show how the wave deforms when it travels away from the source and how the acoustic energy is pumped into higher harmonics from the fundamental frequency.

The Fubini solution is a solution to Equation 1 with no dissipation (δ = 0). Therefore, the difference in amplitudes of the numerical and the analytical solution grows as x comes closer to xsh (compare plots in Figure 3 and Figure 4).

Figure 2: Comparison of the nonlinear numerical solution (blue) with the linear analytical solution (green): the full propagation domain in the top plot and the area around the shock formation distance in the bottom plot. The red vertical line indicates the shock formation distance.

Figure 3: Comparison of model solution (blue) with Fubini nonlinear analytical solution (green) at x = 0.5xsh. The top plot shows the pressure profiles and the bottom plot shows the frequency spectra.

Figure 4: Comparison of model solution (blue) with Fubini nonlinear analytical solution (green) at x = xsh. The top plot shows the pressure profiles and the bottom plot shows the frequency spectra.

Figure 5: Comparison of model solution (blue) with Fay nonlinear analytical solution (green) at x = 3.5xsh. The top plot shows the pressure profiles and the bottom plot shows the frequency spectra.

Notes About the COMSOL Implementation

Mesh and finite element discretization

The mesh is required to resolves the frequency content of the signal. That means resolving higher harmonics along the wave propagation direction. The specified number of harmonics to resolve N0, should therefore contribute to the mesh element size. To accurately resolve the acoustic pressure, use at least second-order (quadratic) elements. Then the mesh element size is defined as dx = c0/(6N0f0).

Time Stepping

The time stepping depends on the maximum frequency to resolve. It is specified in section at the top physics level. In the model it is selected to resolve the desired number of harmonics as fmax = N0f0. The Generalized alpha time stepping method generated as the default transient solver will automatically use an appropriate time step to resolve up to fmax.

Note that when the Nonlinear Acrostics (Westervelt) feature is used the default solver should be regenerated if, for example, a linear model was solved previously. This is to ensure that a proper solve is set up. The addition of the nonlinear feature will trigger the Automatic (Newton) method for solving the boundary value problem.

Shock-capturing stabilization

Since the model discusses a subject of nonlinear wave propagation beyond the shock-formation distance, a special treatment is required to resolve the discontinuities of the shock. The Nonlinear Acoustics (Westervelt) feature contains a built-in Shock-Capturing Stabilization technique available when is enabled in the view. The stabilization is turned off per default as it requires manual tuning.

Whenever the is enabled, an extra nonlinear term

is added to the sound diffusivity, δ. This nonlinear term introduces additional dissipation that is maximal where the acoustic pressure endures discontinuities, that is, at shocks. The parameters κ and q must be tuned so as not to introduce too much or too little dissipation. The choice depends on the material properties and the frequency of the input signal. In this model, κ = 0.01 and q = 1.35. This simple 1D model can be used to tune the stabilization parameters (for other materials and frequencies) for higher dimension models as it is relatively fast to solve.

References

1. M.F. Hamilton and D.T. Blackstock, eds., Nonlinear Acoustics, Academic Press, San Diego, CA, 1998.

2. D.T. Blackstock, Fundamentals of Physical Acoustics, John Wiley & Sons, 2000.

Acoustics_Module/Nonlinear_Acoustics/nonlinear_acoustics_westervelt_1d

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

Load the parameters used in the model from a file. Some of the parameters are presented in Table 1.

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file nonlinear_acoustics_westervelt_1d_parameters.txt.

Now, define the amplitudes Bn used in the Fubini and Fay analytical solutions to Equation 1.

Analytic 1 (an1)

In the toolbar, click

and choose .

In the window for , type Pn_fubini in the text field.

Locate the section. In the text field, type 1/n*besselj(n, n*sigma)*sin(n*omega0*(t - sigma*x_sh/c0)).

In the text field, type sigma, t, n.

Locate the section. In the text field, type 1, s, 1.

In the text field, type 1.

Analytic 2 (an2)

In the toolbar, click

and choose .

In the window for , type Pn_fay in the text field.

Locate the section. In the text field, type 1/sinh(n*(sigma + 1)/Gamma)*sin(n*omega0*(t - sigma*x_sh/c0)).

In the text field, type sigma, t, n.

Locate the section. In the text field, type 1, s, 1.

In the text field, type 1.

Define variables for the linear and nonlinear analytical solutions to the problem. Use the sum() operator with 100 terms to approximate the expression given in Equation 2.

Definitions

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
sigma	x/x_sh		Relative distance
p_fubini	2P0/sigmasum(Pn_fubini(sigma, t, n), n, 1, 100)	Pa	Fubini solution
p_fay	2P0/Gammasum(Pn_fay(sigma, t, n), n, 1, 100)	Pa	Fay solution
p_linear	P0sin(omega0(t - x/c0))	Pa	Linear solution