Jean-François Mercier
email: jean-francois.mercier@ensta.fr

Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

A two-way model for nonlinear acoustic waves in a non-uniform lattice of Helmholtz resonators

Keywords: nonlinear acoustics, solitary waves, Burgers equation, fractional derivatives, diffusive representation, time splitting, shock-capturing schemes

Propagation of high amplitude acoustic pulses is studied in a 1D waveguide connected to a lattice of Helmholtz resonators. An homogenized model has been proposed by Sugimoto (J. Fluid. Mech., 244 (1992)), taking into account both the nonlinear wave propagation and various mechanisms of dissipation. This model is extended here to take into account two important features: resonators of different strengths and back-scattering effects. An energy balance is obtained, and a numerical method is developed. A closer agreement is reached between numerical and experimental results. Numerical experiments are also proposed to highlight the effect of defects and of disorder.

Introduction

We are concerned here with the dynamics of nonlinear waves in lattices. This subject has stimulated researches in a wide range of areas, including the theory of solitons and the dynamics of discrete networks. Numerous studies have been led in electromagnetism and optics [START_REF] Kevrekidis | Nonlinear waves in lattices: past, present, future[END_REF], and many physical phenomena have been highlighted, such as discrete breathers [START_REF] Lazarides | Discrete breathers in nonlinear magnetic metamaterials[END_REF][START_REF] Boechler | Discrete breathers in one-dimensional diatomic granular crystals[END_REF][START_REF] Feng | Discrete breathers in two-dimensional nonlinear lattices[END_REF], chaotic phenomena [START_REF] Grabowski | Wave propagation in a nonlinear periodic medium[END_REF][START_REF] Wan | One-dimensional nonlinear Schrödinger equation: A nonlinear dynamical approach[END_REF], dynamical multistability [START_REF] Wan | Wave transmission in one-dimensional nonlinear lattice: multi stability and noise[END_REF][START_REF] Li | Wave propagation in nonlinear photonic band-gap materials[END_REF] and solitons or solitary waves [START_REF] Li | Scattering properties of solitons in nonlinear disordered chains[END_REF][START_REF] Kartashov | Solitons in Nonlinear Lattices[END_REF]. We focus more specifically on solitary waves, which have been first observed and studied for surface waves in shallow water [START_REF] Russell | Report on Waves, Made to the Meetings of the British Association in 1842-1843[END_REF]. The main feature of these waves is that they can propagate without change of shape and with a velocity depending of their amplitude [START_REF] Remoissenet | Waves Called Solitons: Concepts and Experiments[END_REF]. They have been studied in many physical systems, as in fluid dynamics, optics or plasma physics.

For elastic waves, numerous works have highlighted the existence of solitary waves in microstructured solids [START_REF] Engelbrecht | Nonlinear deformation waves in solids and dispersion[END_REF], periodic structures such as lattices or crystals [START_REF] Chetverikov | Localized nonlinear, soliton-like waves in two-dimensional anharmonic lattices[END_REF][START_REF] Hao | Experiments with acoustic solitons in crystalline solids[END_REF][START_REF] Hess | Solitary surface acoustic waves and bulk solitons in nanosecond and picosecond laser ultrasonics[END_REF], elastic layers [START_REF] Kuznetsov | Soliton-like lamb waves[END_REF][START_REF] Lomonosov | Nonlinear surface acoustic waves: Realization of solitary pulses and fracture[END_REF][START_REF] Mayer | Nonlinear surface acoustic waves: Theory[END_REF], layered structures coated by film of soft material [START_REF] Kovalev | Solitary Rayleigh waves in the presence of surface nonlinearities[END_REF], periodic chains of elastics beads [START_REF] De Billy | On the formation of envelope solitons with tube ended by spherical beads[END_REF][START_REF] Daraio | Tunability of solitary wave properties in one-dimensional strongly nonlinear phononic crystals[END_REF][START_REF] Moleron | Solitary waves in chain of repelling magnets[END_REF] and uniform or inhomogeneous rods [START_REF] Dai | Solitary waves in an inhomogeneous rod composed of a general hyper elastic material[END_REF][START_REF] Khusnutdinova | Fission of longitudinal strain solitary wave in a delaminated bar[END_REF][START_REF] De Billy | Generation of transversal envelope soliton in polymeric and wooden rods[END_REF].

On the contrary, only a few works have dealt with acoustic solitary waves. Such a lack is mainly explained by the fact that the intrinsic dispersion of acoustic equations is too low to compete with the nonlinear effects, preventing from the occurrence of solitons. To observe these waves, geometrical dispersion must be introduced. It has been the object of the works of Sugimoto and his co-authors [START_REF] Sugimoto | Experimental demonstration of generation and propagation of acoustic solitary waves in a air-filled tube[END_REF][START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF][START_REF] Sugimoto | Acoustic solitary waves in a tunnel with an array of Helmholtz resonators[END_REF][START_REF] Sugimoto | Verification of acoustic solitary waves[END_REF], where the propagation of nonlinear waves was considered in a tube connected to an array of Helmholtz resonators. A model incorporating both the nonlinear wave propagation in the tube and the nonlinear oscillations in the resonators has been proposed. Theoretical and experimental investigations have shown the existence of acoustic solitary waves [START_REF] Sugimoto | Experimental demonstration of generation and propagation of acoustic solitary waves in a air-filled tube[END_REF]. We have developed a numerical modeling, and we have successfully compared simulations with experimental data [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF][START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz resonators[END_REF].

One fundamental assumption underlying Sugimoto's model is that all the resonators are the same, which allows to use a homogenization process. The drawback of this approach is that the reflection of an incident wave by a defect (for instance, a resonator different from the others) cannot be considered. Similarly, wave propagation across resonators with variable features cannot be investigated. However, the case of variable resonators is important when studying the influence of manufacturing defects or the influence of aging of the guiding device on the wave propagation.

The aim of this paper is to remedy these limitations, by building a model predicting two-way propagation across variable resonators. For this purpose, we introduce a discrete description of the resonators. A consequence is that the requirement of a long wavelength is no longer necessary. Second, we allow the reflection of waves. Finally, we reformulate the viscothermal losses, which leads to a formulation suitable for an energy balance. The new model writes as a single system coupling three unknowns, the velocity of the right-going wave u + , the velocity of the left-going wave u -and the total excess pressure in the resonators p, induced by both waves:

                   ∂u ∂t + + ∂ ∂x a 0 u + + b (u + ) 2 2 + c a 0 ∂ 1/2 ∂t 1/2 u + -d ∂ 2 u ∂x 2 + = -e(1 -2mp) ∂p ∂t , ∂u ∂t - + ∂ ∂x -a 0 u -+ b (u -) 2 2 + c a 0 ∂ 1/2 ∂t 1/2 u --d ∂ 2 u ∂x 2 - = +e(1 -2mp) ∂p ∂t , ∂ 2 p ∂t 2 + f ∂ 3/2 ∂t 3/2 p + gp -m ∂ 2 (p 2 ) ∂t 2 + n ∂p ∂t ∂p ∂t = h(u + -u -).
(1) The precise meaning of u ± , p and the coefficients in (1) will be detailed along section 2. Some of the coefficients incorporate the individual features of the resonators and vary with x.

The paper is organized as follows. In section 2, the general equations in the tube and in the resonators are given, and the new model ( 1) is derived. In section 3, a first-order formulation is followed: it allows to determine an energy balance and also to build a numerical scheme. In section 4, comparisons with experimental data show that a closer agreement is obtained than with the original Sugimoto's model. Numerical simulations are performed to investigate the properties of the acoustic solitary waves, in particular the robustness to defects and to random disorder. Conclusions and prospects are drawn in section 5. Technical results are given in the appendices.

The new model

In this part, we derive Eq. ( 1). The configuration under study is made up of an air-filled tube connected to uniformly distributed cylindrical Helmholtz resonators (see Fig. 1). The geometrical parameters are the radius of the guide R, the axial spacing between successive resonators D, the radius of the neck r, the length of the neck L, the radius of the cavity r h and the height of the cavity H, which can vary from one resonator to the other. The crosssectional area of the guide is A = π R 2 and that of the neck is B = π r 2 , the volume of each resonator is V = π r 2 h H.

Equations in the tube 2.1.1. General equations

Here we recall briefly the general equations governing the nonlinear acoustic waves in a tube [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF]. The physical parameters are the ratio of specific heats at constant pressure and volume γ; the pressure at equilibrium p 0 ; the density at equilibrium ρ 0 ; the Prandtl number Pr; the kinematic viscosity ν; and the ratio of shear and bulk viscosities µ v /µ. One defines the linear sound speed a 0 , the coefficient of nonlinearity b, the sound diffusivity ν d and the coefficient of attenuation d by

a 0 = γ p 0 ρ 0 , b = γ + 1 2 , ν d = ν 4 3 + µ v µ + γ -1 Pr , d = ν d 2 .
(2) The characteristic angular frequency ω and the parameter of nonlinearity

ε are ω = 2 π a 0 λ , ε = γ + 1 2 u a 0 , ( 3
)
where λ is the typical wavelength imposed by the source, and u is the amplitude of the gas velocity. The following assumptions are made:

• low-frequency regime (ω < ω * = 1.84 a 0 R ), so that only the plane mode propagates [START_REF] Chaigne | Acoustics of Musical Instruments[END_REF];

• weak acoustic nonlinearity: ε 1 [START_REF] Hamilton | Nonlinear Acoustics[END_REF].

Thanks to the first assumption, we can restrict to unidirectional propagation. Under the second assumption of weak nonlinearities, we are allowed to use Riemann invariants along the characteristics, noted with the ± subscript according to the way of propagation. Then, starting from the compressible Navier-Stokes equations, the flow in the tube can be described by the following nonlinear equation for the horizontal velocities u ± [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF]:

∂u ∂t ± + ±a 0 + bu ± ∂u ∂x ± = ± a 0 2ρ 0 F ± + d ∂ 2 u ∂x 2 ± , (4)
with

F ± (x) = 1 A σ(x) (ρv ± n )(x, y(s), z(s))ds. (5)
The diffusivity of sound in the tube is introduced by the coefficient d. F ± is the mean mass flux at position x, with v ± n the small deviation of the component of the velocity normal inward to the inner surface of the tube (v ± n lies in the y -z plane). The integral is taken along the periphery σ(x) of the cross-section at x (see Fig. 2) with s the curvilinear abscissa. F ± is due to the coupling between the main flow u ± and two areas, represented in beige on Fig. 3. It takes into account the connections of the main tube with the resonators, and the wall friction due to the presence of a viscous boundary layer on the boundary of the tube. We note Σ a part of length D (thus connected to only one resonator) of the tube, represented in Fig. 2. We note Σ B the part of Σ connected to the cylindrical resonator, of area B. For simplicity, the junction between a resonator and the main tube is approximated by as a disc. The real shape is more complicated: it is the intersection of two perpendicular tubes of different diameters. It follows [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF]:

ρv ± n = -Q on Σ B , = ρ 0 v ± b on Σ\Σ B , (6)
with

Q = V a 2 0 B (1 -2mp) ∂p ∂t , v ± b = C √ ν ∂ -1/2 ∂t -1/2 ∂u ∂x ± , (7)
where

m = γ -1 2 γ p 0 . (8)
Q is the mass flux density over B into the resonator cavity (see Fig. 2). The reader is referred to the sections 2-1 and 2-2 of [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF] for a detailed derivation of [START_REF] Dai | Solitary waves in an inhomogeneous rod composed of a general hyper elastic material[END_REF]. The expression of Q follows from the conservation laws in the throat and from the adiabatic relation in the cavity, up to the second order.

x z u(x, t)

Q v n v n Figure 3: definition of v ± n .
In beige is represented the viscous boundary layer and the inside of the resonators. v ± n is defined on the boundary of this beige area.

In [START_REF] Dai | Solitary waves in an inhomogeneous rod composed of a general hyper elastic material[END_REF], v ± b is the velocity at the edge of the viscothermal boundary layer. It involves a fractional operator of order -1/2, thus proportional to 1/(i ω) 1/2 in the frequency domain [START_REF] Chester | Resonant oscillations in closed tubes[END_REF]. C is the dissipation in the boundary layer:

C = 1 + γ -1 √ Pr . (9)
These expressions are obtained by solving matching conditions between the boundary layer and the main flow in the guide (section 2-1 of [START_REF] Sugimoto | Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[END_REF]). The Riemann-Liouville fractional integral of order 1/2 of a causal function w(t) (thus w(0) = 0) is defined by

∂ -1/2 ∂t -1/2 w(t) = H(t) √ π t * w = 1 √ π t 0 (t -τ ) -1/2 w(τ ) dτ, ( 10
)
where * is the convolution product in time, and H(t) is the Heaviside step function [START_REF] Podlubny | Fractional Differential Equations[END_REF]. For convenience, Σ is taken as the part x ∈] -r, D -r[, the part Σ B being then centered at x = 0, located in x ∈] -r, r[ (see Fig. 4). Then the flux F ± can be determined at each position x, depending on two cases. Case 1:

x ∈]r, D -r[. Then σ(x) is a disc of radius R and ρv ± n = ρ 0 v ± b on σ. Therefore F ± (x) = 1 A 2π 0 (ρv ± n )(x, R sin θ, R cos θ)Rdθ, = 1 A (2πR) ρ 0 v ± b (x),
where θ is measured with respect to the z axis (as α in Fig. 5(a)).

Case 2:

x ∈] -r, r[. A cross section of the tube at x ∈] -r, r[ is sketched in Fig. 5(a). The integration path σ(x) is made of two parts: (i) a line of length 2 (x) on which ρv ± n = -Q; (ii) a portion of a circle of radius R, between the angles θ = α and θ = 2π -α, on which ρv ± n = ρ 0 v ± b . Using Fig. 5(b) representing Σ B seen from above, it follows that:

(x) = 0 for x ∈]r, D -r[, √ r 2 -x 2 for x ∈] -r, r[. (11)
Using Fig. 5(a), one finds also that α(x) = arcsin ( (x)/R). Then (5) becomes: To sum up, using the values of Q and v ± b from [START_REF] Dai | Solitary waves in an inhomogeneous rod composed of a general hyper elastic material[END_REF], one deduces that for all x ∈] -r, D -r[:

F ± (x) = 1 A 2π-α(x) α(x) (ρ 0 v ± b )(x, R sin θ, R cos θ)Rdθ - (x) -(x) Q(x)ds , = 2 A [π -α(x)] Rρ 0 v ± b (x) -Q(x) (x) . (a) (b) z y
↵(x) R `(x) x z r ` (x)
a 0 2ρ 0 F ± = c(x) ∂ -1/2 ∂t -1/2 ∂u ∂x ± e(x)(1 -2mp) ∂p ∂t , (12)
where we have introduced the functions

   c(x) = c 0 1 - 1 π arcsin (x) R , e(x) = e 0 (x), (13)
with the coefficients

c 0 = C a 0 √ ν R , e 0 = V ρ 0 a 0 AB , (14)
and V the volume of the resonator at x = 0. This result has been obtained for a unit part of the tube located in x ∈ ] -r, D -r[. Now, we consider the entire tube, in which the resonators are

D H n 1 H n H n+1 x n+1 x n 1
x n 2r Figure 6: side view of the guide connected to Helmholtz resonators. The i-th cell is highlighted in grey.

centered at 6). We note V i the volume (and H i the height) of the i-th resonator at

x i = iD, i = 1, 2, • • • (see Fig.
x i . In the i-th cell x ∈ [x i -D/2, x i +D/2]
(represented in grey on Fig. 6) and based on (12), Eq. ( 4) takes the form

∂u ∂t ± + ∂ ∂x ±a 0 u ± + b (u ± ) 2 2 ∓c(x) ∂ -1/2 ∂t -1/2 ∂u ∂x ± -d ∂ 2 u ∂x 2 ± = ∓e(x)(1-2mp) ∂p ∂t . (15)
The coefficients c(x) and e(x) are defined by

   c(x) = c 0 1 - 1 π arcsin i (x) R , e(x) = e i i (x), (16)
where

i (x) =    0 for r < |x -x i | ≤ D 2 , r 2 -(x -x i ) 2 for |x -x i | ≤ r, (17)
c 0 is defined in [START_REF] Fouque | Wave propagation and time reversal in randomly layered media[END_REF], and

e i = V i ρ 0 a 0 AB . (18)

Reformulation of the viscothermal losses

The fractional integrals in [START_REF] Grabowski | Wave propagation in a nonlinear periodic medium[END_REF] have two drawbacks: they prevent from obtaining an energy balance, and the dispersion relation is singular at zero frequency [START_REF] Lombard | Diffusive approximation of a time-fractional Burgers equation in nonlinear acoustics[END_REF]. It is therefore preferable to replace them by fractional derivatives. Based on the linear approximation (∂u ± /∂t) ± a 0 (∂u ± /∂x) ≈ 0, one has

∂ -1/2 ∂t -1/2 ∂u ± ∂x ≈ ∓ 1 a 0 ∂ 1/2 ∂t 1/2 u ± . ( 19
)
The rigorous derivation of this step is detailed in Appendix A. Injecting [START_REF] Kartashov | Solitons in Nonlinear Lattices[END_REF] in [START_REF] Grabowski | Wave propagation in a nonlinear periodic medium[END_REF] yields finally the first two equations of (1).

Equation in the resonators

The excess pressure p in a resonator of volume V satisfies the equation [START_REF] Monkewitz | The response of Helmholtz resonators to external excitation. Part 1. Single resonators[END_REF][START_REF] Monkewitz | The response of Helmholtz resonators to external excitation. Part 2. Arrays of slit resonators[END_REF]:

∂ 2 p ∂t 2 + f ∂ 3/2 p ∂t 3/2 + gp -m ∂ 2 p 2 ∂t 2 + n ∂p ∂t ∂p ∂t = gp , (20)
where

f = 2 √ ν r L L e , g = a 2 0 B L e V , n = V B L e ρ 0 a 2 0 , (21)
and where p is the excess pressure at the mouth of the tube. The coefficients m in [START_REF] Daraio | Tunability of solitary wave properties in one-dimensional strongly nonlinear phononic crystals[END_REF] and n in (21) describe nonlinear processes in the resonators. The semi-empirical coefficient n accounts for the jet loss resulting from the difference in inflow and outflow patterns [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF][START_REF] Sugimoto | Verification of acoustic solitary waves[END_REF]. These nonlinear processes have to be included to get a good agreement with the experimental measurements [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF]. The Caputo fractional derivative of order 3/2 is obtained by applying (10) to ∂ 2 p/∂t 2 . Corrected lengths have been introduced: L = L + 2 r accounts for the viscous end corrections, and the corrected length L e = L + η accounts for the end corrections at both ends of the neck, where η ≈ 0.82 r is determined experimentally [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF].

In the one-way model [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF], a linear approximation is used to link the right-going fields in the tube: p + = (γp 0 /a 0 )u + . In a symmetric way, p -= -(γp 0 /a 0 )u -. Here p is induced by both the right-going waves and the left-going waves; assuming linearity gives

p = p + + p -= γp 0 a 0 (u + -u -). ( 22
)
Injecting ( 22) in ( 20) leads to

∂ 2 p ∂t 2 + f ∂ 3/2 p ∂t 3/2 + gp -m ∂ 2 p 2 ∂t 2 + n ∂p ∂t ∂p ∂t = h(u + -u -), (23)
with h = g(γ p 0 /a 0 ). For resonators of variable volume, Eq. ( 23) becomes, for |x -

x i | ≤ r (i = 1, 2, • • •) ∂ 2 p ∂t 2 + f ∂ 3/2 p ∂t 3/2 + g(x)p -m ∂ 2 p 2 ∂t 2 + n(x) ∂p ∂t ∂p ∂t = h(x)(u + -u -), (24)
with

g(x) = a 2 0 B L e V i , n(x) = V i B L e ρ 0 a 2 0 , h(x) = γ p 0 a 0 g(x). (25)
One recognizes the third equation in [START_REF] De Billy | On the formation of envelope solitons with tube ended by spherical beads[END_REF]. The functions g(x), h(x) and n(x) depend on x because the resonators have different heights, but they are constant on each resonator. They just depend on i, the index of the resonators. The equation ( 24) is not solved outside the resonators (that is on r < |x -

x i | ≤ D/2) contrary to Sugimoto's model recalled in (B.6b).
In Appendix B, we show that the original model of Sugimoto [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF] can be recovered from the new model [START_REF] De Billy | On the formation of envelope solitons with tube ended by spherical beads[END_REF]. Note also that (1) can be easily extended to the case of resonators not periodically placed; only the values of x i need to be changed.

First-order system

In this part, Eq. ( 1) is formulated as a first-order system, which enables us to determine an energy balance and to build a numerical scheme.

Diffusive approximation

The fractional derivatives in (1) are nonlocal in time and they rely on the full history of the solution, which is numerically memory-consuming. An alternative approach is based on a diffusive representation of fractional derivatives, and then on its approximation. This method has already been presented in [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF][START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz resonators[END_REF] and we just recall it briefly. The half-order integral [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives[END_REF] of w(t) is written

∂ -1/2 ∂t -1/2 w(t) N =1 µ ϕ (t), (26)
where the diffusive variables ϕ (t) = ϕ(t, θ ) satisfy the ODE

   ∂ϕ ∂t = -θ 2 ϕ + 2 π w, ϕ (0) = 0. ( 27
)
The approximation [START_REF] Li | Scattering properties of solitons in nonlinear disordered chains[END_REF] follows from the approximation of an integral thanks to a quadrature formula on N points, with weights µ and nodes θ , which are issued from an optimization process.

To get fractional derivatives of orders 1/2 and 3/2, we differentiate (26) in terms of t. Then we deduce:

∂ 1/2 w ∂t 1/2 = ∂ ∂t ∂ -1/2 w ∂t -1/2 N =1 µ ∂ϕ ∂t (t) = µ -θ 2 ϕ + 2 π w . (28)
Similarly, the derivative of order 3/2 is written

∂ 3/2 w ∂t 3/2 = ∂ ∂t ∂ 1/2 w ∂t 1/2 N =1 µ ∂ξ ∂t (t), (29)
where we have introduced ξ = ∂ϕ /∂t. The diffusive variable ξ (t) satisfies the following ODE (derivative of ( 27)):

   ∂ξ ∂t = -θ 2 ξ + 2 π ∂w ∂t , ξ (0) = 0. ( 30
)
The initial condition is obtained using w(0) = 0, since we consider causal data. Thanks to these diffusive approximations, Eq. ( 1) can be written as the following system

                                               ∂u ∂t ± + ∂ ∂x ±a 0 u ± + b (u ± ) 2 2 = - c a 0 N =1 µ -θ 2 ϕ ± + 2 π u ± + d ∂ 2 u ∂x 2 ± ∓ e(1 -2mp)q, ( 31a
)
∂p ∂t = q, ( 31b
)
∂q ∂t = h(u + -u -) -gp -f N =1 µ -θ 2 ξ + 2 π q + m ∂ 2 p 2 ∂t 2 -n |q| q, ( 31c
)
∂ϕ ± ∂t = -θ 2 ϕ ± + 2 π u ± , = 1 • • • N, (31d)
∂ξ ∂t = -θ 2 ξ + 2 π q, = 1 • • • N. (31e)
Contrary to [START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz resonators[END_REF], all the initial conditions are null: u ± = 0, p = 0, ∂p ∂t = 0. If this were not the case, then the diffusive representation of the 3/2 derivative would imply non-null initial conditions on ξ in [START_REF] Lombard | Diffusive approximation of a time-fractional Burgers equation in nonlinear acoustics[END_REF]. The interested reader is referred to [START_REF] Lombard | Diffusive approximation of a time-fractional Burgers equation in nonlinear acoustics[END_REF] for additional details on this topic.

Energy balance

The system (31) is suitable to define an energy and to prove the energy decreasing, in an infinite medium and for smooth solutions:

Result 1. Let E = E + 1 + E - 1 + E 2 and K = K + 1 + K - 1 + K 2 , (32)
with

                                         E ± 1 = 1 2 R (u ± ) 2 + π 2 c a 0 N =1 µ θ 2 (ϕ ± ) 2 dx, (33a)
E 2 = 1 2 R eg h p 2 + e h (1 -2 m p) q 2 + π 2 ef h N =1 µ θ 2 ξ 2 dx, (33b)
K ± 1 = π 2 R c a 0 N =1 µ ∂ϕ ± ∂t 2 dx + R d ∂u ± ∂x 2 dx, (33c)
K 2 = R e n h q 2 |q| - m n q + π 2 e f h N =1 µ ∂ξ ∂t 2 dx . ( 33d
)
We neglect the 2mp terms in the tube, ie in the evolution equations of u ± in (1). Then the following energy balance holds:

dE dt = -K. ( 34
)
The proof is reported in Appendix C. Four remarks are raised by [START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF].

• E ± 1 and K ± 1 involve quantities in the tube, notably the kinetic energy (u ± ) 2 ; E 2 and K 2 involve quantities in the resonators, notably the potential energy with terms proportional to p 2 and q 2 .

• We have not succeeded in obtaining a proof by accounting for the 2mp term in the advection equations. Nevertheless, the hypothesis 2mp 0 is consistent with Sugimoto's work, where only the influence of m in the resonators equation has been considered: see (B.6a). It is also consistent with the hypothesis of a weak nonlinear regime. Indeed, as shown in [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF][START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz resonators[END_REF]:

2mp ≈ (γ -1) u + -u - a 0 , (35)
which is lower than 1 under the hypothesis of weak nonlinearity |u ± | a 0 .

• The term E in (32) is positive if µ > 0 and 1 -2mp > 0. The first condition is imposed when the coefficients of the diffusive representation are determined during the optimisation process [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF]. The second condition is satisfied in the weak nonlinear regime.

• The term K in ( 32) is positive if µ > 0 and m < n. The first condition has been already discussed. The second condition reads m/n = γ-1

2 BLe V < 1,
where BL e is the resonator neck volume, and V is the volume of the resonators. For the experimental configuration under study, this ratio is lower than 1 and thus the condition is satisfied.

Numerical scheme

In this part, we present the numerical resolution of (31) with null initial conditions. Source terms s ± (t) at x s ± model the acoustic sources of rightgoing waves and left-going waves: u + (x s + , t) = s + (t), u -(x s -, t) = s -(t). For numerical purpose, it is necessary to write the evolution equations as a firstorder system in time. To do so, the term ∂ 2 p 2 /∂t 2 in Eq. ( 31c) is expanded and leads to

(1 -2mp) ∂q ∂t = h(u + -u -) -g p -f N =1 µ -θ 2 ξ + 2 π q + 2mq 2 -n |q| q.
The (4 + 3 N ) unknowns for the counter-propagating waves are gathered in the vector

U = u + , u -, p, q, ϕ + 1 , • • • , ϕ + N , ϕ - 1 , • • • , ϕ - N , ξ 1 , • • • , ξ N T . ( 36
)
Then the nonlinear systems [START_REF] Lomonosov | Nonlinear surface acoustic waves: Realization of solitary pulses and fracture[END_REF] is written in the form

∂ ∂t U + ∂ ∂x F(U) = G ∂ 2 ∂x 2 U + S(U), (37)
where F is the flux function

F = a 0 u + + b (u + ) 2 2 , -a 0 u -+ b (u -) 2 2 , 0, 0, • • • , 0 T , (38)
and where the source term S is

S =                               - c a 0 N =1 µ -θ 2 ϕ + + 2 π u + -e(1 -2mp) q - c a 0 N =1 µ -θ 2 ϕ -+ 2 π u -+ e(1 -2mp) q q 1 1 -2mp h(u + -u -) -g p -f N =1 µ -θ 2 ξ + 2 π q + 2mq 2 -n |q| q -θ 2 ϕ + + 2 π u + , = 1 • • • N -θ 2 ϕ -+ 2 π u -, = 1 • • • N -θ 2 ξ + 2 π q, = 1 • • • N                               . (39)
As soon as m = 0 and n = 0, S(U) is no longer a linear operator (m = 0 = n has been considered in [START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz resonators[END_REF]). The Jacobian matrix ∂F ∂U in ( 38) is diagonalizable with real eigenvalues: a 0 + b u + , -a 0 + b u -and 0 with multiplicity 3 N + 2, which ensures propagation with finite velocity. These eigenvalues do not depend on the quadrature coefficients µ and θ . The diagonal matrix G = diag(d, d, 0, • • • , 0) incorporates the volume attenuation.

Let us describe briefly the numerical procedure; additional details can be found in [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF][START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz resonators[END_REF]. To compute the N quadrature coefficients µ and θ in [START_REF] Remoissenet | Waves Called Solitons: Concepts and Experiments[END_REF], we use a nonlinear optimization with the positivity constraints µ ≥ 0 and θ ≥ 0 [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF]. In order to integrate the system (37), a grid is introduced, with a uniform spatial mesh size ∆x and a variable time step ∆t n . The computations are done with N x = 1000 grid nodes. The approximation of the exact solution U(x j = j ∆ x, t n = t n-1 + ∆t n ) is denoted by U n j . A stability analysis [START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz resonators[END_REF] leads to the restriction on the time step

β = a (n) max ∆t n ∆x 1 + 1 Pe ≤ 1. ( 40
)
In the CFL condition [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF], Pe = a

a (n) max = a 0 + b max j [max(|u n+ j |, |u n- j |)]. (41)
In practice, we choose β = 0.95. To treat the system (37), a Strang splitting is used [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction[END_REF], ensuring both simplicity and efficiency: the original equation ( 31) is split in a propagative equation and a forcing equation, which are solved successively with adequate time increments. The propagative part is solved by a standard second-order TVD scheme (a finite-volume scheme with flux limiters) for nonlinear hyperbolic PDE [START_REF] Leveque | Numerical methods for conservation laws[END_REF] combined with a centered finitedifference approximation [START_REF] Lombard | Numerical modeling of nonlinear acoustic waves in a tube connected with Helmholtz resonators[END_REF].

Results

γ p 0 (Pa) ρ 0 (kg/m 3 ) P r ν (m 2 /s) µ v /µ 1.

Test 1: comparison with experiments

The first result is based on the experimental measurements detailed in [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF]. The reader is referred to this paper for a complete description of the experimental and numerical setups. The physical and geometrical parameters are given in table 1. The explosion of a balloon at x = 0 generates a shock wave, which propagates along a tube of length 6.15 m. This wave interacts with resonators regularly spaced from x = 0.2 m up to the end of the tube. Two models are compared: the one-way model with constant coefficients (ie the original Sugimoto's model used in [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF]) and the two-way model with variable coefficients developed here. The computations are done on 5000 grid nodes. A preliminary numerical study has been performed to verify that this discretization is sufficiently fine to capture the converged solution, at the scale of the figures. The same remark holds for all the forthcoming tests. Fig. 7 (a-b) compares the excess pressure p , defined in Eq. ( 22), simulated by the models or measured experimentally. Two heights are considered: For both models, a good agreement is observed between the simulated data and the experimental results. The best agreement is obtained with the new model, especially for larger resonators heights. It is observed that higher resonators yield smaller waves with a larger support. Moreover, the arrival time of the central peak is around 10 ms in (a) and 14 ms in (b), which indicates that higher resonators generate slower waves, as predicted by the theory [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF].

Test 2: amplitude dependence of the velocity

From now on, we investigate the properties of the two-way model thanks to numerical experiments. As a second test, we consider the spatial domain [0, 10] m, discretized on 5000 grid nodes with N = 4 memory variables. A time-variable punctual source at x s = 5 m is used to generate right-going and left-going waves:

u ± (x s , t) = A ± G(t), (42)
where G is a causal Gaussian pulse with G(0) = 0 and with a central frequency f c = 650 Hz. The amplitudes are A + = 20 m/s and A -= A + /2 = 10 m/s. Figure 8(a,b) display snapshots of p , defined in Eq. ( 22), at different times, for two resonators heights: H = 0 cm (a), and H = 2 cm (b). In the first case without resonators, the viscothermal losses are insufficient to prevent the occurence of shocks in finite time, as analysed theoretically in [START_REF] Sugimoto | Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[END_REF]. On the contrary, resonators create smooth waves (b). From Figure 8, one estimates the velocity V of the nonlinear waves. Figure 9 shows the results obtained by varying the amplitude A + of the forcing, from 10 m/s to 100 m/s, and for two heights H of the resonators. It is observed that V increases linearly with A + , which is a typical signature of solitary waves. At a given forcing, V is also greater for smaller height H, as in the first test of section 4.1.

Test 3: interaction with defects

In this test, the source of right-going waves is localized at the left edge of the computational domain (x = -5 m). The amplitude of the forcing is A + = 20 m/s; from now on, A -= 0 m/s. The height of the resonators is H = 2 cm, except for the resonator at x = 0 m, where H = 0.1 cm. Figure 10 displays snapshot of p at different times. The position of the defect is denoted by a vertical dotted line. At t = 10 ms, the incident wave has not yet crossed the defect (figure 10(a)) and exhibits a smooth shape already seen in the second test. As a fourth and last test, we examine the effect of a random disorder on the propagation of nonlinear waves. The domain [-15, 15] m is discretized on 15000 grid nodes, ensuring 10 grid nodes inside each resonator. The random height of each resonator is uniformly distributed between H min and H max , which gives rise to random values of the coefficients e, g, n and h in [START_REF] De Billy | On the formation of envelope solitons with tube ended by spherical beads[END_REF]. Figure 11 illustrates two configurations. The first configuration has a mean height H = 2 cm, with H min = 0.1 cm and H max = 3.9 cm, and a forcing amplitude of right-going waves A + = 20 m/s (a,c). The second configuration has a mean height H = 7 cm, with H min = 1 cm and H max = 13 cm, and a forcing amplitude A + = 100 m/s (b,d). A zoom between 0 and 1 m illustrates the values of H in 11 resonators (in blue), denoted by vertical solid lines (a,b).

Figure 11(c,d) shows a snapshot of p at t = 80 ms without disorder (blue line) and with a random selection of H (red line). With randomness, a coda is observed behind the main wavefront. The amplitude and the location of the peak are only slightly modified, which indicates the robustness of the nonlinear wave when perturbed by some noise. It constitutes an additional signature of solitary waves.

Conclusion

The goal of this paper was to build a model describing the propagation of nonlinear wave across a variable lattice of Helmholtz resonators. Contrary to previous works [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF][START_REF] Sugimoto | Verification of acoustic solitary waves[END_REF][START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF], the new model takes into account the individual features of each resonator and the backscattering effects. Doing so has enabled us to study the effect of disorder or of an individual defect on the propagation of nonlinear waves.

The numerical experiments have highlighted various properties of the new model. First, the Helmholtz resonators prevent from the occurence of shocks and yield large-amplitude smooth structures. Moreover, the velocity of waves increases linearly with the amplitude, which confirms the terminology of "acoustic solitary waves". These properties were already included in the original Sugimoto's model [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF], but here a better agreement is obtained with experimental results. Second, diffraction effects have been displayed when a local variation of the height of the resonator is introduced in the lattice. Third and last, robustness to disorder has been observed by performing simulations in random lattices.

This work motivates new experimental investigations. The setup used in [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF] enables to examine whether the effects seen in figures 10-11 are confirmed experimentally. New theoretical and numerical [START_REF] Xiu | Numerical methods for stochastic computations: a spectral method approach[END_REF] developments to study the propagation in a stochastic medium are also motivated, for instance concerning the robustness of acoustic solitary waves in random media [START_REF] Fouque | Wave propagation and time reversal in randomly layered media[END_REF]. Lastly, adequate choices of defects in the lattice could produce localized modes, generalizing the results presented in the linear case in [START_REF] Sugimoto | Localized mode of sound in a waveguide with Helmholtz resonators[END_REF].

Acknowledgement

The authors thank O. Richoux who gave access to his experimental measurements. Agnès Maurel is also greatly acknowledged for her support. We thank also the anonymous Reviewers for their constructive remarks.

Appendix A. Transformation of the fractional integral

Here we prove [START_REF] Kartashov | Solitons in Nonlinear Lattices[END_REF]. To simplify, we restrict to a right-way propagation in [START_REF] Grabowski | Wave propagation in a nonlinear periodic medium[END_REF], and we note u = u + and p = p + :

∂u ∂t + ∂ ∂x a 0 u + b u 2 2 -c ∂ -1/2 ∂t -1/2 ∂u ∂x -d ∂ 2 u ∂x 2 = -e(1 -2mp) ∂p ∂t . (A.1)
with the coefficients defined in ( 8) and ( 16). The first step is to write (A.1) in a dimensionless form. Given a central frequency f 0 , we define the characteristic wavelength λ = a 0 /f 0 . We introduce the non-dimensional quantities, indexed by a tilde: x = λx, t = t/f 0 , u = u 0 ũ and p = p 0 p, where u 0 is the characteristic velocity, that will be determined later. We get:

∂ ũ ∂ t + ∂ ∂ x ũ + b ũ2 2 - c ∂ -1/2 ∂ t-1/2 ∂ ũ ∂ x - d ∂ 2 ũ ∂ x2 = -ẽ(1 -2 mp) ∂ p ∂ t , (A.2) with b = M b = M γ + 1 2 , c = c a 0 λ a 0 , d = d a 0 λ , ẽ = ep 0 u 0 = V 2ADγM , m = mp 0 = γ -1 2γ . (A.3) M = u 0 /c
= ẽ ⇔ M = V ADγ(γ + 1) . (A.4)
To find an approximation of (A.2), we introduce the small parameter ε = max( b, c, d, ẽ). Then we define new quantities, defined with hats, by the relation v = εv, with v = b, c, d or e, such that all the quantities with a hat are at most equal to one. Then, starting from (A.2), we get:

∂ ũ ∂ x = - ∂ ũ ∂ t -ε bũ ∂ ũ ∂ x + ε d∂ 2 ũ ∂ x2 -εê(1 -2 mp) ∂ p ∂ t + εĉ ∂ -1/2 ∂ t-1/2 ∂ ũ ∂ x , = - ∂ ũ ∂ t -ε bũ ∂ ũ ∂ x + ε d∂ 2 ũ ∂ x2 -εê(1 -2 mp) ∂ p ∂ t + εĉ ∂ -1/2 ∂ t-1/2 - ∂ ũ ∂ t -ε bũ ∂ ũ ∂ x + ε d∂ 2 ũ ∂ x2 -εê(1 -2 mp) ∂ p ∂ t + εĉ ∂ -1/2 ∂ t-1/2 ∂ ũ ∂ x , = - ∂ ũ ∂ t -ε bũ ∂ ũ ∂ x + ε d∂ 2 ũ ∂ x2 -εê(1 -2 mp) ∂ p ∂ t -εĉ ∂ 1/2 ũ ∂ t1/2 + O(ε 2 ),
where we have used

∂ 1/2 ũ ∂ t1/2 = ∂ -1/2 ∂ t-1/2 ∂ ũ ∂ t .
Therefore (A.2) can be approximated by:

∂ ũ ∂ t + ∂ ∂ x ũ + b ũ2 2 - c ∂ 1/2 ũ ∂ t1/2 - d∂ 2 ũ ∂ x2 = -ẽ(1 -2 mp) ∂ p ∂ t , (A.5)
where neglected terms are small, of order ε 2 . Using the values of Table 1, M = 0.21 is obtained from (A.4), in agreement with the experiments made in [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF], for f 0 = 500 Hz. From this Mach value, we deduce finally b = 0.22 = ẽ, c = 0.010, d = 8.2410 -8 , m = 0.14. (A.6)

Coming back to real quantities yields the first two equations in [START_REF] De Billy | On the formation of envelope solitons with tube ended by spherical beads[END_REF].

• in the tube, the coefficient of adiabatic nonlinearity 2mp is neglected;

• an averaged description of the geometry is used: for identical resonators of volume V , a continuous approximation of the tube geometry is introduced and a mean value of the flux is used:

F ± = 1 D D x=0 F ± (x)dx. (B.1)
Using ( 13) and ( 14 Also, thanks to the approximation for R r:

arcsin √ r 2 -x 2 R √ r 2 -x 2 R . (B.3)
Using ( 13) and ( 14 Using these two ingredients, the system (1) degenerates in two families of uncoupled equations: which corresponds to the model proposed in [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF][START_REF] Sugimoto | Verification of acoustic solitary waves[END_REF] and studied in [START_REF] Richoux | Generation of acoustic solitary waves in a lattice of Helmholtz resonators[END_REF]. In (C.1c), the derivative of p 2 is modified, using

           ∂u ∂t ± + ∂ ∂x ±a 0 u ± + b (u ± ) 2 2 ∓ c ∂ -1/2 ∂t -1/2
µ ∂ϕ ± ∂t + d ∂ 2 u ± ∂x 2 ∓
∂ 2 ∂t 2 (p 2 ) q = ∂ ∂t (p q 2 ) + q 3 . (C.3)
This relation is easy to check by expanding both sides of the equality. Next we multiply (C.1c) by q, and using (C.1b), (C.2) and (C.3) leads to

q ∂q ∂t = h(u + -u -) q -g p q -f N =1 µ ∂ξ ∂t q + m ∂ 2 p 2 ∂t 2 q -n |q| q 2 , = h(u + -u -) q -g p ∂p ∂t - π 2 f N =1 µ ∂ξ ∂t + θ 2 ξ ∂ξ ∂t + m ∂ ∂t (p q 2 )
-nq 2 |q| -m n q .

(C.4) Therefore, we get

(u + -u -)q = 1 2 ∂ ∂t g h p 2 + 1 h (1 -2 m p) q 2 + π 2 f h N =1 µ θ 2 ξ 2 + π 2 f h N =1 µ ∂ξ ∂t 2 + n h q 2 |q| - m n q .
(C.5)

Besides, (C.1a) is multiplied by u ± and integrated in space. After summation and integration by parts (the data are compactly supported), we get

R u + ∂u + ∂t + u -∂u - ∂t dx = - R c a 0 N =1 µ u + ∂ϕ + ∂t + u -∂ϕ - ∂t dx - R d ∂u + ∂x 2 + ∂u - ∂x 2 dx - R e (u + -u -) q dx.
(C.6) Thanks to the relations (C.2) and (C.5), the previous equation is simplified and the conclusion follows.

Figure 1 :

 1 Figure 1: sketch of the guide connected with Helmholtz resonators.

Figure 2 :

 2 Figure 2: perspective view of the tube. The unit portion Σ of the tube is connected to one resonator through the surface Σ B . σ is a path along the boundary of a cross section of the tube.

2. 1 . 2 .Figure 4 :

 124 Figure 4: portion of length D of the tube. View from the side (a). View from above (b).

Figure 5 :

 5 Figure 5: (a) σ(x) for x ∈] -r, r[; (b) view from above of the resonators in x ∈] -r, r[.

  ∆x/2 d ≈ 10 5 is the discrete Péclet number, with the maximal velocity

Figure 7 :

 7 Figure 7: Test1. Time history of the excess pressure p at x r = 2.8 m. Comparison between the models and the experiments. (a): H = 7 cm; (b): H = 13 cm.

Figure 8 :

 8 Figure 8: Test 2. Snapshots of the excess pressure p generated by a source point at x = 5 m, at various instants. (a): H = 0 cm (no resonators); (b): H = 2 cm.

Figure 9 :

 9 Figure 9: Test 2. Velocity of the nonlinear waves versus the forcing amplitude A + for different heights of the resonators.

Figure 10 :

 10 Figure 10: Test 3. Snaphsots of p at t = 10 ms (a) and t = 25 ms (b). All the resonators have the height H = 2 cm, except at x = 0 (vertical dotted line), where H = 0.1 cm.

Figure 10 (

 10 Figure 10(b) shows p at t = 25 ms, after the interaction with the defect. The amplitude of the transmitted wave (in red) has been slightly modified, compared to the perfect case without defect (in blue). More important, a left-going wave has been generated by the defect. This effect cannot be

4. 4 .

 4 Test 4: propagation in a random medium

Figure 11 :

 11 Figure 11: Test 4. Left: H = 2 cm (a-c); right: H = 7 cm (b-d). Top: discretizations of the random heights. The grid nodes are denoted by blue points; the vertical dotted lines denote the positions of the resonators. Bottom: snapshots of p at t = 80 ms in the unperturbed medium (red) and in the random medium (blue).

  0 is the characteristic Mach number. The typical value of M is obtained by equalizing b and ẽ: b

  eq, (C.1a)∂p ∂t = q, (C.1b) ∂q ∂t = h(u + -u -) -gp -f N =1 µ ∂ξ ∂t + m ∂ 2 p 2 ∂t 2 -n |q| q, (C.1c) ∂ϕ ± ∂t = -θ 2 ϕ ± + 2 π u ± , = 1, • • • , N, (C.1d)

Table 1 :

 1 Physical parameters of the air at 15 • C, and geometrical parameters of the tube with resonators.

	403	10 5	1.177	0.708 1.57 10 -5 0.60
	R (m) D (m) r (m)	L (m) r h (m)	H (m)
	0.025	0.1	0.01	0.02	0.0215	0.02, 0.07 or 0.13

Appendix B. Comparison with Sugimoto's model

Here we prove that the original model of Sugimoto [START_REF] Sugimoto | Propagation of nonlinear acoustic waves in a tunnel with an array of Helmholtz resonators[END_REF] can be recovered from the new model [START_REF] De Billy | On the formation of envelope solitons with tube ended by spherical beads[END_REF]. In this aim, we introduce three ingredients:

• a restriction to only one-way propagations: for right-going waves (u -= 0), we note p + the associated pressure p. Similarly, for left-going waves (u + = 0), we note p -the pressure p;

Appendix C. Energy balance

Here we prove the result 1. The first step is to rewrite the system (31), eliminating the term 2m|p| 1 in (31a) and using [START_REF] Lombard | Diffusive approximation of a time-fractional Burgers equation in nonlinear acoustics[END_REF], to get the following system: