Numerical investigation of acoustic solitons

Acoustic solitons can be obtained by considering the propagation of large amplitude sound waves across a set of Helmholtz resonators. The model proposed by Sugimoto and his coauthors has been validated experimentally in previous works. Here we examine some of its theoretical properties: low-frequency regime, balance of energy, stability. We propose also numerical experiments illustrating typical features of solitary waves.


Introduction
Solitons are nonlinear waves with large amplitude and constant profile, resulting from the competition between nonlinearity and dispersion.They occur in many physical area, such as fluid mechanics (Kortewegde Vries equations), electromagnetism and optics (Klein-Gordon equations) [1].In acoustics, the intrinsic dispersion is too low compared to the nonlinearity to produce solitons.Thus additional geometric dispersion must be considered to observe acoustic solitons.It was the basis of a series of works of Sugimoto and coauthors [6,7], where the propagation of shock waves was investigated in a tube connected to an array of Helmholtz resonators.A mathematical model was proposed, as well as a theoretical analysis and a comparison with experimental data.
Sugimoto's work was extended in two means.In [3], a time-domain numerical model was proposed to incorporate efficiently the fractional derivatives modeling linear viscothermic losses.In [5], comparisons with experimental results were proposed.It was shown that nonlinear attenuation in the resonators had also to be incorporated for describing accurately the experiments.
The goal of the present contribution is to analyse further the full Sugimoto's model with fractional derivatives and nonlinear attenuation, recalled in section 2.1.In the low-frequency regime, corresponding to the experimental conditions, the evolution equations tend to a Korteweg-de Vries equation with an additional nonlinear term.Therefore one expects solitons nonlinearly attenuated.Both for mathematical and numerical purposes, the fractional model is transformed by means of a diffusive representation (section 3.1.).Doing so allows to analyse the energy balance and the stability of the model (sections 3.2.and 3.3.).Lastly, two sets of numerical experiments are proposed in section 4., showing that the waves have the typical features of solitary waves.

Sugimoto's equations
The configuration is depicted in figure 2 of [3].The wavelengths are much larger than the distance between two resonators, so that the latter are described by a continuous distribution.One-dimensional propagation is assumed.The variables are the velocity of gas u and the excess pressure in the resonators p. Considering the right-going wave, one writes [6] The PDE (1a) describes the nonlinear wave propagation (coefficients a and b) in the tube.The losses in the tube are introduced by c (viscothermic losses at the wall) and by d (volume attenuation).The ODE (1b) describes the oscillations in the resonators.In the latter, the losses are introduced by f (viscothermic losses), by m and by n (nonlinear attenuation due to turbulence).Coupling between (1a) and (1b) is ensured by e and h.See [5] for the expression of all these coefficients.A fractional integral of order 1/2 (c) and a fractional derivative of order 3/2 ( f ) are introduced.These non-local operators are tackled with in section 3.1.

Low-frequency approximation
Under the hypothesis of weak nonlinearity, ∂ u/∂ x in (1a) is replaced by −(1/a) ∂ u/∂t in the terms with coefficients b, c and d.The resulting system is written in the (T, X ) coordinates, where T is a nondimensional retarded time, X is a non-dimensional slow space variable: where u max is the magnitude of the gas velocity at the initial time, ω is a characteristic wave frequency, and γ is the ratio of specific heats at constant pressure and volume.Introducing the reduced variables U and where p 0 is the pressure at equilibrium, one obtains the system This system generalizes the equations (2-5) and (2-6) of [7] to the case of nonlinear losses (terms with M and N).As shown in [6], β is negligible compared to δ r and δ R .We consider waves with characteristic frequencies much smaller than the natural frequency of the resonators ω e , so that Ω = (ω e /ω) ≫ 1.In this case, the dispersion analysis performed in [3] indicates that the viscothermal losses are small.Moreover, the volume of the resonators is large compared to that of the necks, so that M ≪ N is neglected.Consequently, the low-frequency regime Ω ≫ 1 yields the simplified system From (5b), one obtains Injecting ( 6) in (5a) gives: Neglecting the second-order terms in 1/Ω and introducing the new unknown When nonlinear attenuation in the resonators is neglected (N = 0), equation ( 7) recovers the Korteweg-de Vries equation (2-35) of [6], which allows the propagation of solitons.Solitons are also expected to exist for small N values, but with a decrease of amplitude.

Evolution equations
A diffusive approximation of the non-local in time fractional operators in (1) is followed here [4].The half-order integral of a function w(t) can be written where the diffusive variable φ satisfies the local-in-time ordinary differential equation In (8), φ (t, θ ℓ ) = φ ℓ (t); µ ℓ and θ ℓ are the weights and nodes of the quadrature formula.Their computation is detailed in [5].A similar derivation is applied to the 3/2 derivative in (1), involving the diffusive variable ξ .Injecting these diffusive approximations in (1) yields the following system of evolution equations The (3 + 2 N) unknowns are gathered in the vector Then the nonlinear system (10) can be written in the form

Energy balance
Based on the system (10), we define Assuming smooth solutions (no shock) and c = 0, then one obtains It follows two remarks.First, if the weights of the diffusive approximation are positive µ ℓ > 0, then E is a quadratic definite positive form, which thus defines an energy.In practice, we determine these weights by an optimization procedure with constraint of positivity [5].Second, if the coefficient of nonlinear attenuation satisfies m = 0, then dE /dt < 0: the energy decreases, and the model is well-posed.In practice, this hypothesis is reasonable, since m ≪ n.
Let us finally examine the assumed hypotheses.With shocks, the wave motion is irreversible and additional terms of dissipation must be accounted for in (14).On the other hand, the hypothesis c = 0 is not physical but required for technical purpose: up to now, we have not found an energy if c = 0.

Stability analysis
The system (12) is solved by a splitting technique [5]: one successively solves the PDE with adequate time steps.Here we examine the stability of both stages.First, ∂ F ∂ U has real eigenvalues {a + b u, 0 2 N+2 } and is diagonalizable.Consequently, ( 15) is hyperbolic when G = 0.In practice, G introduces a parabolic regularization, and the problem remains well-posed.
From cases (ii) and (iii), it follows that the spectral radius of the Jacobian satisfies ρ(T) = |λ N | > θ 2 N ≫ 1.As a consequence, (16) must be solved by an implicit scheme [5].We conjecture that this stiffness of T still holds for nonlinear attenuation (m = 0 or n = 0) and for the general case (10).This justifies the splitting strategy.

Numerical results
In this part, we examine whether the solution of the Sugiomoto's model (1) has the typical features of solitons.In test 1, one investigates the dependence of the velocity upon the amplitude of the forcing.In test 2, we simulate the interaction between two waves.The physical and geometrical parameters are given in [5].Two values of the resonators height are considered: H = 2 cm and H = 7 cm.This parameter influences the resonance angular frequency of the Helmholtz resonators (ω e in section 2.2.) and the parameters K and N in (7).The waves are generated by imposing the value of the velocity in (10) at x = 0.A Gaussian with amplitude A is chosen for this purpose: ( The central frequency is f 0 = 1/t 0 = 650 Hz.The standard deviation τ is chosen so that u(0, 0) = u(0, 2t 0 ) = A/1000.A set of 10 receivers is distributed uniformly on the computational domain.Seismograms are built from the time signals stored.The positions of the maximal value of u at each receiver is detected and allows to estimate the celerity V of the wave.An example for H = 7 cm and A = 100 m/s is given in figure 1-(a).After a transient regime (offsets 0 and 1), a smooth structure emerges despite the nonsmoothness of the evolution equations (10).The amplitude of the wave decreases along propagation, due to the loss mechanisms.Lastly, small amplitude waves are observed before the main wave front.The same procedure is followed by varying A from 10 m/s to 100 m/s.The evolution of V in terms of A is illustrated in figure 1-(b).The linear increase of V with A is clearly observed.Greater values of V are obtained for smaller value of H.These two observations confirm the theoretical analysis performed in [7].

Interaction of two solitons
A Gaussian pulse with small amplitude followed by a taller pulse are generated.Due to its higher amplitude, the latter travels faster, allowing an interaction between the two soliton waves.Figure 2 presents the results of this experiment.In the inviscid case (a), we observe that the two waves interact in a manner analogous to classical solitons [2]: after the waves separate, each one has again the form of a solitary wave, though shifted in location from where they would be without interaction (denoted by crosses).When the attenuation mechanisms are accounted for (b), a similar observation can be done, even if the observation is not so clear due to the smoothing of waves.

Conclusion
In this contribution, we have studied some properties of the full Sugimoto's model with nonlinear attenuation.Theoretical analysis has shown that the coefficient m can produce problems (increase of energy, loss of stability).Since this coefficient has a negligible practical influence, we propose to remove it.Numerical experiments have allowed to examine situations difficult to reproduce experimentally.They have shown that typical features of solitons are maintained despite the nonlinear attenuation mechanisms.

4. 1 .Fig. 1 .
Fig. 1.Test 1. (a): example of seismogram.The vertical dotted lines represent the location of the maximum at each receiver.The inclined red line denotes the trajectory of these maxima; its slope yields the velocity of the wave.(b): velocity of the waves in terms of the forcing amplitude.