OSCILLATING NONLINEAR ACOUSTIC SHOCK WAVES

. We investigate oscillating shock waves in a tube using a higher order weakly nonlinear acoustic model. The model includes thermoviscous ef- fects and is non isentropic. The oscillating shock waves are generated at one end of the tube by a sinusoidal driver. Numerical simulations show that at res- onance a stationary state arise consisting of multiple oscillating shock waves. Oﬀ resonance driving leads to a nearly linear oscillating ground state but su- perimposed by bursts of a fast oscillating shock wave. Based on a travelling wave ansatz for the ﬂuid velocity potential with an added 2’nd order poly- nomial in the space and time variables, we ﬁnd analytical approximations to the observed single shock waves in an inﬁnitely long tube. Using perturbation theory for the driven acoustic system approximative analytical solutions for the oﬀ resonant case are determined.


1.
Introduction. In cylindrical tubes nonlinear standing acoustic waves may appear as oscillating shock waves. If acoustical resonators are driven at or near one of their resonance frequencies, the amplitude of the acoustical field inside the resonator becomes very high. Within the theory of lossless linear acoustics it is well-known that the amplitude goes to infinity at the resonance frequencies. However, in real situations the amplitude at resonance becomes finite due to dissipative and nonlinear effects. In typical situations dissipation is small and the amplitude encountered will be large. Accordingly, nonlinear effects come into play and may give rise to oscillating shocks inside the resonator. Much of the basic theory for nonlinear standing waves is due to Chester (1964) [2], who studied resonant oscillations in closed tubes. A study of nonlinear standing waves based on a perturbation approach applied to the Kuznetsov equation is given in the book by Enflo and Hedberg (2002) [1].
It is a well-known fact that nonlinear acoustical waves may generate a flow field in which the particle velocities are not simply sinusoidal, and a pattern of timeindependent vortical flows or steady circulations is often found in the body of compressible media. Such flow patterns are known as acoustic streaming (Nyborg, 1998) [15] and (Hagsäter et al., 2007(Hagsäter et al., , 2008 [5,6].
Here we investigate nonlinear standing waves in a one-dimensional acoustical resonator having one end closed, and one end periodically driven. The analysis is based on numerical simulations of a higher order nonlinear acoustical model previously presented in [16,17]. In the literature many different weakly nonlinear models have been proposed and investigated focusing on various specific physical properties of fluids. Some prominent examples are found in [12,3,4,9,11]. Solutions of these nonlinear acoustic models are normally not easy to determine analytically, but travelling waves have been found by Jordan in [8] and furthermore soliton like solutions of the nonlinear acoustic Lagrangian-averaged Euler-α model has been determined in [11]. Some further mathematical treatments can be found in the review of Kaltenbacher [10].
If we drive the oscillator at one of its resonance frequencies, and choose the amplitude of the periodic excitation and the thermoviscous dissipation parameter appropriately, the solution consists of a number of standing shocks inside the oscillator. We demonstrate that driving an acoustic wave at low frequencies far from resonance results in a ground state standing wave superimposed an oscillating shock wave. The appearance of the shock wave depends on the phase of the driver. In the half period where the driver compresses the fluid an oscillating shock wave travels forth and back in the tube and in the half period where the driver expands the fluid no shock wave appears. Parallel to numerical simulations results we have used a new travelling wave ansatz for shock waves superimposed second degree polynomials in the space and time variables. The ansatz is for the fluid velocity potential. For the off resonance driving case we present a solution ansatz and perturbation method for a slowly varying driving force. This approach is compared to a Galerkin expansion approach using four eigenmodes.

2.
A higher order nonlinear acoustic wave equation. A number of simplified model equations for weakly nonlinear acoustic wave propagation have been derived in the literature [16,17,3,4,9,11]. Here we use a model for acoustic waves in a Newtonian, viscous and heat conducting gas. Invoking the dynamical equations for the fluid motion, continuity, the heat transfer of temperature end entropy together with an equation of state, one can derive the one dimensional plane wave model [16] whereΨ =Ψ(X, T ) is the fluid velocity potential at position X and at time T . The fluid flow velocity U (X, T ) is given by the potential through U = −Ψ X . Subscripts X and T on dependent variables asΨ denote partial derivatives with respect to the space and time variables, respectively. The fluid pressure P is given by P =Ψ T , and the speed of sound is c 0 in the undisturbed fluid. The parameter γ is the specific heat under constant pressure divided by the specific heat under constant volume. The pressure P of the fluid is expanded to second order in the fluid density and first order in entropy from the equilibrium state. In this expansion γ is proportional to the second order term in the pressure [13]. The parameterb is the diffusivity of sound and this term depends on shear viscosity, bulk viscosity and the heat conductivity coefficient [7,16]. In other wordsb is a loss term. The equation in (1) can be written in dimensionless variables by introducing the scaling The independent space variable x is dimensionless and the constant k x carries the dimension meter (m). The constant k t has the dimension seconds (s) and t is dimensionless. Finally, k ψ has the dimension m 2 /s leaving ψ(x, t) dimensionless. By introducing the dimensionless damping parameter b and choosing we obtain the dimensionless model equation Note that k t can be chosen freely. If the numerical part of k t without the unit s is selected identical to the numerical part of 1/c 2 0 (without unit) the Eq. (4) becomes identical to Eq. (1) for c 0 = 1m/s. If b = 0 in Eq. (4) the remaining parts of the model constitute a Hamiltonian system with the Lagrangian From this Lagrangian it is easy to derive the associated Hamiltonian, and the Euler Lagrange variational equation leads to Eq. (4) with b = 0. In the following we consider acoustic wave propagation in a cylindrical tube of normalized length ℓ assuming plane waves. The x-axis is along the tube center axis. At the left end (x = 0) of the tube a sound generator is mounted and at the right end (x = ℓ) we have a fixed wall. The corresponding boundary conditions become u(0, t) = −ψ x (0, t) = D sin(ωt) and u(ℓ, t) = −ψ x (ℓ, t) = 0 .
Here D is the driver amplitude and ω is the driving frequency. The initial conditions corresponds to a fluid at rest, that is ψ(x, 0) = 0 and ψ t (x, 0) = 0. In the simulations below we solve Eq. (4) using these initial conditions and solve numerically until steady state has emerged. We observe that Eq. (4) do not possess time reversal symmetry. That is changing the sign of t do not leave the model equation Eq. (4) unchanged. We shall later see that this lack of time reversal symmetry has striking consequences for the solutions. In order to solve Eq. (4) numerically we have transformed the governing equations into a set of two first order coupled differential equations In the above system we have introduced u 1 (x, t) = ψ(x, t) and u 2 (x, t) = ψ t (x, t). The corresponding boundary conditions at x = 0 are u 1,x (0, t) = −D sin(ωt) and u 2,x (0, t) = −Dω cos(ωt). At the right hand side of the tube we have u 1,x (ℓ, t) = 0 and u 2,x (ℓ, t) = 0. The system has been discretized in space using second order finite differences keeping time as a continuous variable. The resulting system of ordinary differential equations has been integrated numerically using the Matlab Runge Kutta 4-5 algorithm. This method is often referred to as the method of lines.  (4), ω = 2π corresponds to the eigenfrequency of the second harmonic. This means we drive the nonlinear equation (4) at a resonance frequency. However, due to damping and the nonlinear terms the emerging steady state solution is two oscillating shock waves, travelling forth and back in opposite directions. Generally, exciting the n'th eigenmode results in n standing shocks. Our numerical results agree qualitatively with the results presented in Chester (1964) [2] and Enflo and Hedberg (2002) [1].
The right panel of Fig. 1 shows a countour plot of the oscillating shock waves. Even though they are highly nonlinear we do not observe a phase shift during collision of the shocks as for solitons. The numerical results exhibit shock waves propagating with velocities equal to unity. We have calculated the average pressure B(t) defined by as function of time t. We have found a slow linear pressure build up as time progress from B(0) = 0 to B(100) = 0.04. We believe this pressure build up is a result of acoustic streaming caused by the nonlinearities in the model. Here the fluid cannot stream out at x = ℓ and hence the reaction is a slight pressure build up.

3.2.
Driving at nonresonant frequency. In Fig. 2 we show the fluid velocity u(x, t) as function x and t obtained numerically for a nonresonant driving frequency at ω = 0.1 and with the driving amplitude D = −0.125. The other parameter values used are the same as in Fig. 1, i.e. cavity length ℓ = 1, damping paramter b = 5 · 10 −4 and nonlinearity parameter γ = 1.4. We find the results striking. The slowly varying ground state oscillation is superimposed by a faster oscillating shock wave. The shock wave propagates at the velocity one and a magnification of the shock wave oscillation is shown in the right panel of Fig. 2. The simulations reveal that the shock wave oscillations appear for decreasing u(0, t) corresponding to compression of the fluid and disappears for increasing u(0, t) corresponding to decompression.

4.
Analysis of a sloping shock wave ansatz. The first travelling wave solution determined analytically for a higher order nonlinear acoustics model we believe is the one presented by Jordan in [8] for the Kuznetsov equation [12]. We shall follow the ideas in [8] with modifications adopted to the standing shocks observed in the numerical simulations presented in Fig. 1. These waves resemble sums of traveling shocks and functions that are nearly linear in x and t, i.e. traveling shocks with sloping lines on each side of the shock. Guided by this fact we search for exact solutions by introducing the following ansatz [17] ψ(x, t) = Ψ(ξ) + a 1 x + a 2 t + a 3 xt + a 4 where ξ = x − vt. We note that Eq. (8) is an extension of the solution ansatz used in reference [16] for thermoviscous shock solutions, where a 3 = a 4 = a 5 = 0. In the following analysis we shall find exact solutions, which are not the same shock waves as observed in the numerical simulations presented in the previous section. We actually find new exact shock wave solutions with propagation speeds less than unity in contrast to the shock waves generated above travelling with velocities equal to unity. Note that the ansatz (8) includes second order terms in x and t. Taking the derivatives of Eq. (8) with respect to x and t, respectively, yields where Φ = −Ψ ′ and prime denotes differentiation with respect to ξ. Inserting Eq. (8) into Eq. (4) yields where Equation (10) depends explicitly on both x and t. However, by choosing α 1 = α 2 = α 3 = α 4 =0 we eliminate this dependence. Solving Eq. (11) with vanishing α i , i = 1, 2, 3, 4 gives where a 5 can be chosen arbitrarily. By inserting Eqs. (12) into Eq. (10), and invoking Φ = −Ψ ′ we obtain By integrating Eq. (14) once with respect to ξ, we get It is convenient to change variables according to In terms of these variables Eq. (15) takes the form where σ = sign(a 5 ) and δ = 2 |a 5 |b 3v 2 − 1 .
Thus, Eq. (17) describes the motion of a mechanical particle according to in the potential The parameter ζ plays the role of time but in contrast to common mechanical systems where time is always positive, t ∈ (0, ∞), we have here ζ ∈ (−∞, ∞). Therefore when δ > 0 the particle loses energy for ζ > 0 (positive damping) and gains energy for ζ < 0 (negative damping). While for δ < 0 the particle loses energy for ζ < 0 (positive damping) and gains energy for ζ > 0 (negative damping). It is interesting to note that Eq. (19) can be written in an equivalent conservative form. To this end we introduce instead of the function χ(ζ) a new functionχ(ζ): By inserting Eq.(21) into Eq. (19), we get Thus Eq. (22) describes the motion of a mechanical particle according to in a "time-dependent" potential The Lagrangian for Eqs. (23) and (24) is given by Let us first consider the case when δ = 0 corresponding to v 2 = 1/3 in Eq. (18). Note that this case occurs in media with γ = 1, see Eq. (12). Now the equations (17) and (22) become identical and have the form The solution of Eq. (26) can be obtained in implicit form where the constant E has the meaning of energy of the effective particle. We shall consider the cases a 5 > 0 (σ = 1) and a 5 < 0 (σ = −1) separately.

4.2.
The case σ = −1. The potential function V (χ) for the case σ = −1 is presented in Fig. 3. The effective particle oscillates in the potential well. The value E = 0 separates two regimes. For E > 0 the particle hits the vertical wall χ = 0 (note that in accordance with definition we have χ ≥ 0). When E < 0 the function never vanishes, that is χ = 0. Taking into account the definition (16), one can conclude that the only physically interesting case is for E < 0, because for E > 0 and where W (k, x), k = 0, ±1, ±2, . . ., is the Lambert-W function at branch point k, that is the solution of the equation W e W = x [14]. The period of the oscillations is given by As it follows from Eq. (30) while for E → 0 the period slowly diverges as (See Fig. 4). For 0 < E + 1 4 ≪ 1 one can approximate the potential function V (χ) as a series By inserting Eq. (33) into Eq. (27), one can get where h j (j = 1, 2, 3), h 3 > h 2 > h 1 are roots of the cubic equation  Fig.6. Both in Fig. 6 and Fig.  5 the function − dΨ dζ is plotted for ζ ∈ (− T 2 , T 2 ). Note that γ = 1 corresponds to isothermal acoustic propagation.

4.3.
The case σ = 1. In this case the potential function has the form presented in Fig. 7. In the frame of the mechanical analogy this means that for χ > 1 the particle moves downhill and its coordinate rises indefinitely. In the interval χ ∈ (0, 1) the particle oscillates in the well with a wall at the vertical line χ = 0.  Fig. 8 obtained by using the Matlab solver ode45. In Fig. 8(a) we show the solution for a 5 = −0.05. It is seen that the behavior presented in Fig. 6 qualitatively agree with the results obtained by solving Eq. (13) numerically shown in Fig. 8(a). Figure 8(b) shows another solution with a 5 = 0.1.
Once Φ is computed numerically we insert the data into Eqs. (9) in order to obtain u(x, t) and η(x, t). Next we determine the function Ψ by numerically integrating Φ with respect to ξ and by insertion into Eq. (8) we obtain ψ(x, t). In order to verify  5. Analysis of the non resonant case. In the following section we shall analyse the non resonant case in more detail of the problem with an acoustic driver at the  Fig. 8(a). The solution is shown at five instants during t ∈ [0, 10].  Figure 10. Continued from Fig. 9. The initial and boundary conditions are obtained from the solution given in Fig. 8(b).
left hand side of a tube. In order to do so we introduce the ansatz where the function ϕ(x, t) satisfies the following equation with homogeneous boundary conditions Here and in what follows prime denotes derivatives with respect to t. Let us consider a slowly varying driving force at the boundary x = 0 In this case one can neglect the terms with U ′ and U ′′ in Eq. (37) and the equation takes the form Let us present the function ϕ(x, t) as a sum where φ has a zero mean value and thus By inserting Eq. (41) into Eq. (40), we get where and It is seen that Eq. (46) has the solution φ(x, t) = 0. In this case G = 0 and Eq.
(44) has a solution where W (0, x) is the principal branch of the Lambert-W function. The coefficient A is determined by the initial value of Φ ′ (t) The mean pressure in the tube is given by the expression for ℓ = 1. As an example let us consider the case of periodic driving In this case The mean pressure is a periodic function of time. The amplitude of the oscillations depends on the material parameter δ, ω and the intensity of the driving D. In the full numerical simulation of the nonresonant case of the model in Eq. (4) we observed a ground state wave superimposed an oscillating shock wave, see Fig. 2.
The above perturbation analysis has not captured the shock wave. The same holds for the Galerkin approach to 4'th order in Eq. (53). The analysis in this section is valid alone for the ground state oscillation.

6.
Conclusion. Numerical studies of acoustical driven resonators revealed the presence of oscillating shock waves inside the resonator similar to nonlinear kink anti-kink propagation in long Josephson junctions modelled by the perturbed sine-Gordon equation. Driving at the frequency corresponding to the n'th linear excitation mode in the linearized model of (4) leads to n oscillating fully nonlinear shocks. However, we believe there is an upper limit in n corresponding to the number of shock waves, which can fit into the cavity. This is also observed in long Josephson junctions.
For the nonresonant driving case full numerical simulations revealed the surprising result, that the excited nearly linear ground state is superimposed oscillating shock waves. These shock waves governed by Eq. (4) travel with a velocity equal to unity and arise in bursts. For decreasing u(0, t) at the left boundary of the resonator oscillating shocks appear and they disappear for increasing u(0, t). We believe the oscillating shocks are connected to the lack of time reversal symmetry of Eq. (4), however, this assertion needs further investigation.
The particular wave pattern observed in the acoustical resonators, lead us to derive new exact generalized shock wave solutions using a localized travelling wave ansatz for the fluid velocity potential superimposed by a 2'nd degree polynomial in the space and the time variables. However, these exact solutions travel at a speed of v = 1/ 3(2 − γ), which is less than unity observed for the travelling shocks in the driven resonator. The reason is that the analytical solutions are for an infinite cavity.
Finally, we have developed a perturbation analysis for the nonresonant driving case invoking a periodic factor on a 2'nd order polynomial in the space variable matching the correct driving boundary condition at the left hand side of the acoustic resonator. Approximate analytical expressions have been obtained, which agree well with an approach based on a 4'th order Galerkin method. However, the perturbation method cannot capture the oscillating shock waves observed in the full numerical simulations of the weakly nonlinear accoustic model in (4).