Propagation of Quasi-plane Nonlinear Waves in Tubes

This paper deals with possibilities of using the generalized Burgers equation and the KZK equation to describe nonlinear waves in circular ducts. A new method for calculating of diffraction effects taking into account boundary layer effects is described. The results of numerical solutions of the model equations are compared. Finally, the limits of validity of the used model equations are discussed with respect to boundary conditions and the radius of the circular duct. The limits of applicability of the KZK equation and the GBE equation for describing nonlinear waves in tubes are discussed.


Introduction
Propagation of nonlinear sound waves in waveguides is a very interesting physical problem.In the case that nonlinear waves travel in a gas-filled waveguide, we can observe phenomena such as nonlinear distortion, nonlinear absorption, diffraction, lateral dispersion, boundary layer effects, etc.All these phenomena can be described by means of the complete system of the equations of hydrodynamics, see e.g.[1]- [4]: the Navier-Stokes momentum equation, the continuity (mass conservation) equation, the heat transfer (entropy) equation and the state equations.Unfortunately, there is no known general solution of this system of equations, and numerical solutions bring many problems regarding stability of the solutions and their time consumption.Consequently, it is sensible to simplify the fundamental system of equations if we ignore some phenomena or consider some of them weak.This simplification leads to the derivation of model equations of nonlinear acoustics.
There have been a number of papers devoted to various aspects of the propagation of nonlinear waves in waveguides.The viscous and thermal dissipative effects on the nonlinear propagation of plane waves in hard-walled ducts are treated for instance in papers [5], [6], [7].The authors deal with the dependence of the frequency on the dissipative and dispersive effects induced by the acoustic boundary layer.Experimental results focused on propagation of finite-amplitude plane waves in circular ducts are presented in [8] and [9].Here, a very good agreement is demonstrated between experimental data and results obtained by means of the Rudnick decay model for the fundamental harmonic.Burns in [10] obtained a fourth-order perturbation solution for finite-amplitude waves.However, his expansion breaks down for large times because it contains secular terms.He took into account dissipation, but he neglected mainstream dissipation with respect to boundary dissipation.Keller and Millmann in [11] found the solution of the model equation for inviscid isentropic fluids, where they used perturbation expansion adapted to eliminate secular terms and determined the nonlinear wavenumber shift for dispersive modes.In [12], Keller utilized the results from [11].He rewrote the results in a form that is useful near the cutoff frequency, in order to show that the cutoff frequencies and resonant frequencies of modes in acoustic waveguides of finite length depend upon the mode amplitude.Nayfeh along with Tsai in [13] presented the nonlinear effects of gas motion as well as the non-linear acoustic lining material properties on wave propagation and attenuation in circular ducts.They obtained a second-order uniformly valid expansion by using the method of multiple scales.These authors presented [14], where they investigated nonlinear propagation in a rectangular duct whose side walls were also acoustically treated by means of the method of multiple scales.Ginsberg in [15] dealt with nonlinear propagation in rectangular ducts.He determined by an asymptotic method the nonlinear two-dimensional acoustic waves that occur within a rectangular duct of semi-infinite length as the result of periodic excitation.In [16], he utilized the perturbation method of renormalization to study the nonlinearity effect on a hard-walled rectangular waveguide.Nonlinear wave interaction in a rectangular duct was also investigated by Hamilton and TenCate in [17] and [18].Multiharmonic excitation of a hard-walled circular duct was treated by Nayfeh in [19].He used the method of multiple scales to derive a nonlinear Schrödinger equation for temporal and spatial modulation of the amplitudes and the phases of waves propagating in a hard-walled duct.
Foda presented his work in [20], which is concerned with the nonlinear interactions and propagation of two primary waves in higher order modes of a circular duct, each at an arbitrary different frequency and finite amplitude.He used the renormalization method to annihilate the secular terms in the obtained expression.If we take no diffraction effects into account, we can use the generalized Burgers equation to describe nonlinear plane waves in circular ducts.It is the Burgers equation that is supplemented by the term that represents boundary layer effects, see [21]- [25].The generalized Burgers equation enables the description of dissipative and dispersive effects that are caused by the boundary layer.Asymptotic and numerical solutions of this equation were presented by Sugimoto [26].In [27], there is a description of the solution of the Burgers equation in the time domain.The approximate solutions of the generalized Burgers equation for harmonic excitation in both the preshock and the postshock region are presented in [28].In the case that only weak diffraction effects are considered, the KZK equation can be used, see [22], [25], [29].However, the boundary layer is incorporated into the boundary condition in [22] and [25].

Basic equations
From the equations of continuum mechanics we can obtain in quadratic approximation an equation describing weakly nonlinear waves in thermo-viscous fluids ( ) ( Here j is the velocity potential, t is time, c 0 is the small-signal sound speed, b is the dissipation coefficient of the medium, r 0 is the ambient density, g is the adiabatic index, equal to the ratio of the specific heat at constant pressure c p to that at constant volume, c V , and L is the second-order Lagrangian density, which is given as ( ) The symbol D represents the Laplacian operator and the symbol Ñ is the gradient operator.The operators are defined for axisymmetric waves in the cylindrical coordinates by where z is the coordinate along the axis of the tube, r is the coordinate perpendicular to the axis of the tube, e z is the unit vector along the z axis and e r is the unit vector along the r axis.We shall assume that the profile of the wave distorts slowly as it propagates in the positive z direction.The term "slowly" means that the wave must travel a distance of many wavelengths l for its profile to be distorted significantly.Then it is reasonable to seek a solution having this functional form ( ) where t is retarded time and m is the small dimensionless parameter which is given by the acoustic Mach number -the ratio of the particle velocity amplitude v m to the small-signal sound speed c 0 .If we use the new coordinate system (t, z 1 , r 1 ), then transformation of the partial derivatives yields Substitution of the partial derivatives (6) into the operator (4) gives Here ^denotes directions transverse to propagation.We can write the linear wave equation directly from Eq. ( 1) If we use the operator (10) in Eq. ( 11) for the divergence operator, then after neglecting a quadratic order term we obtain If we suppose the transverse partial velocity then we can omit the last term in Eq. ( 12): Equation ( 14) represents the plane progressive wave impedance relation.
When conditions ( 13) and ( 14) are satisfied we can suppose ).This enables us to modify Eq. ( 1) Substituting relation (14) to the second-order term of Eq. ( 15) we get where After using the transformation of derivatives ( 6) we can neglect the cubic and higher terms.This yields the well- . The KZK eqatuion was derived on the condition that v z ~m and v ^~m 2 (i.e.v z v ^).v z and v ^are connected by the irrotationality condition of the velocity field If we do not take into account the wall friction we have to use for the solution of the KZK the following boundary conditions where R 0 is the tube radius and ( ) F r, t represents the prescribed sound field at z = 0.
If we assume the wall friction in the case that nonlinear waves propagate in a rigid tube, then a thin boundary layer appears near the tube wall.Within the boundary layer, the velocity component in the direction of the tube axis decreases from a mainstream value to zero at the tube wall.The boundary layer affects the acoustic waves not just near the walls but in the entire volume.The boundary layer causes both energy dissipation and wave dispersion.
If the boundary layer is assumed to have a small displacement effect on the mainstream, the following relation can be obtained, see [24] Here B is the coefficient which is given as where n is the kinematic viscosity, Pr is the Prandtl number, the fractional derivative of the order1 2 in Eq. ( 22) represents the following integrodifferential operator ( ) , , and b T is the thermal expansion coefficient of the fluid b r ¶r where T is temperature.If we suppose a perfect gas then b g Relation ( 22) is valid on condition that where d is the boundary layer thickness where w is an angular frequency.
Differentiating Eq.( 14) with respect to r, we get the equation (the irrotationality condition) If we take into account Eqs. ( 29) and ( 22), we can modify boundary condition (21) for tubes with wall friction as follows , , d ¶t After space-averaging over a whole tube cross section, if we take into account boundary condition (30), the KZK equation is reduced to the generalized Burgers equation where ¢ = B B R 0 .The solution of the linearized GBE for the boundary condition ( ) ( ) where a is the attenuation coefficient for the classical thermo--viscous loss mechanism and the attenuation coefficient a b represents the losses due to the wall friction The ratio of the two attenuation coefficients can be written as It is obvious that d l >1for high frequencies but R 0 l 1 because the first of the conditions (27) has to be satisfied.This means that the classical thermo-viscous losses dominate for high frequencies in comparison to the boundary layer losses.
Consequently, the condition d l 1 cannot be satisfied for higher wave form harmonics.

Conditions for the quasi-plane wave in a tube
The condition for validity of the KZK equation in a tube is given by the relation between the values of the longitudinal and transverse velocity components v v ^.This condition follows from the derivation of the KZK equation in the tube given above, see section Basic equations.
The transverse velocity component relates to transverse diffraction.The value of v ^depends on wave frequency, the non-quasi planar boundary condition, the tube radius and the boundary layer.
Supposing that the tube is narrow, the influence of the boundary layer is significant.Let us assume that the source frequency is lower than the cut-off frequency, thus only plane--wave mode propagation is possible.Then the diffraction is given by the influence of the boundary layer only.When we now gradually increase the frequency, the transverse component of velocity v ^also gradually rises.As soon as the frequency exceeds the cut-off frequency of the tube, the transverse mode arises, thus the effect of transverse diffraction increases considerably and the transverse component of the velocity significantly rises.Oscillations appear on the propagation curves for the harmonic amplitudes which were obtained by means of a numerical solution of the KZK equation and consequently the KZK equation loses its accuracy.This means that in the case of a narrow tube we can study wave propagation only below the cut-off frequency.In addition, it is not possible to study wave propagation if we use a boundary condition of any distinct velocity distribution along the tube radius (non-quasi planar boundary condition).
The influence of the boundary layer is less significant in wide tubes.The amplitudes of transverse components of velocity are much lower than in the case of narrow tubes.Therefore the propagation can be studied far above the cut-off frequency even with distinct velocity distribution along the tube radius for the boundary condition.
However, the KZK equation cannot be used for the boundary condition with distinct velocity distribution along the radius when the source frequency is less than the cut-off frequency for a given tube.This restriction holds for any tube radius.
If the used frequency is less than the cut-off frequency, then all higher modes are evanescent and only the plane wave propagates along the waveguide.
However, the KZK equation cannot describe this process because the transverse component of the velocity is large and we cannot consider the KZK equation to be valid.
The cut-off frequency of the first symmetric mode in a circular waveguide is given by the formula where k 1 is the radial wave number of the first mode and 38317 is the first zero point of the Bessel function derivative: ( ) . Now, let us suppose that the time dependency of the dimensionless longitudinal component near the source is given as where x = r R 0 and q wt = is the dimensionless retarded time, which means that we suppose a mono-frequency harmonic source.Here ( ) g 0 x represents the function describing the distribution of the longitudinal component along the radius.Then for the transverse component of the velocity we have ( ) where W is the dimensionless transverse component of the velocity normalized in the same way as the longitudinal component V.
Assuming wave propagation with frequency ¢ f equal to one half of the cut-off frequency we can see that W is equal to Taking into account the fact that the maximum of the derivative of the commonly used velocity distributions (Gauss, Fermi) at the boundary of a tube is of the order O(1), we can see that the values of v ^and v z are comparable.Therefore the validity condition of the KZK equation is not fulfilled.We can say that the KZK equation is valid only for a moderate distribution along the tube radius below the cut-off frequency (for instance, a plane wave which is affected by the boundary layer).
We often read (e.g., in [4], [29]) the condition kR 0 1 from which follows the relation between wave vector components k k ^.This means that a sound wave is quasi-plane and diffraction effects are weak in tubes.With regard to the text above we can say that the condition kR 0 1 represents only a sufficient, but not a necessary, condition for the use of the KZK equation as a model equation for nonlinear waves propagating in tubes.

Numerical solution
Equation ( 17) was solved in the frequency domain.When the solution is periodic in time (i.e., the sound source is periodic in time), it can be expressed in the form of Fourier series.For a numerical solution it was necessary to truncate the infinite number of terms in the Fourier series to * terms.Then the solution was sought in the form where s b w = z v c zm 0 2 and x = r R 0 are dimensionless coordinates.The finite number of terms in series (41) causes instability in the numerical solution.When the solution is approaching the region of shock wave formation, all harmonics are excited and the energy flow stops with the last harmonic, i.e., the M-th harmonic.This effect causes the Gibbs oscillations in the numerical solution.These oscillations were damped by the method described, e.g., by Fenlon [31].Each harmonic was multiplied by the coefficient Y n given by ( ) where H is the frequency damping coefficient.This causes the additional artificial attenuation of the solution.The value of H was chosen so that Gibbs oscillations practically did not arise.During the assembling of the numerical schema we proceeded from the well known Bergen code [32], which was extended by terms including the boundary layer influence.Simple iteration by the LU decomposition method was used to calculate the solution in the next layer.If we use LU decomposition, then the calculations take approximately three times longer than in the case of calculation by means of simple iteration.However, LU decomposition is not restricted by the condition on the step size in the direction of propagation.Application of this method is very convenient in the case when it is necessary to carry out a very large number of iterations by means of the simple iteration method.For example, we need 500000 steps for the simple iteration method and only 400 steps for LU decomposition for a tube 4 mm in radius.
The Burgers equation was solved in the frequency domain by means of the standard Runge-Kutta method of the fourth order.

A narrow tube
The presented distributions represent the deformation of a plane wave under the influence of the boundary layer for a narrow tube, where the boundary layer effect is very important.The calculations are done for particle velocity amplitude v m = 2.76 m×s -1 , tube radius R 0 = 0.004 m, and for the following six frequencies 1 kHz, 10 kHz, 40 kHz, 60 kHz, 100 kHz and 500 kHz.The cut-off frequency of the first symmetric mode for the given tube is approximately equal to 52 kHz.The first 100 harmonics were used in the numerical calculation, the numerical step in the direction of the wave propagation was Ds = 0.01 for the Burgers equation.The value of the numerical step for the KZK equation was determined from this step, so the total number of steps was identical for the GBE and KZK equation.The KZK equation was solved with the plane boundary condition.The numerical step in the direction of the tube radius was chosen Dx = 0.05.The calculations are realized for the case of propagation in air.The source oscillated harmonically with one frequency.For all calculations we used the numerical attenuation coefficient H = 40.
Figures 1-3 contain three curves: curve 1 is the numerical solution of the GBE, curve 2 is the longitudinal component of the obtained by the numerical solution of the KZK equation, and curve 3 represents the transverse velocity component which is calculated from the longitudinal velocity across the tube radius using eq.( 18).The values of curve 3 are multiplied 100 times and calibrated in the same way as the longitudinal component of the velocity to enable the two velocity components to be compared easily.Curves 2 and 3 are depicted for the radius x = 05 . .Owing to the fact that the resultant space evolutions correspond to practically plane waves, the choice of this radius is not important.Figures 1-3 contain the first harmonic.
The graphs illustrate the speculation about the validity of the KZK equation for a narrow tube, where the influence of the boundary layer is dominant.We see that with increasing frequency the transversal velocity component progressively rises.We can observe the significant growth of the transversal velocity component when the cut-off frequency is exceeded.The first oscillations of the transversal velocity component are noticeable in the case of 60 kHz.These oscillations still grow and they are conspicuous when the frequency is equal to 100 kHz.The longitudinal component of the velocity was scarely affected, as can be seen from the comparison with the velocity from GBE.The further increase in frequency causes the lon- gitudinal component of velocity also to become disrupted, as can be seen for frequency equal to 500 kHz.These variations are readily noticeable, especially for the first harmonic, and they are weaker for higher harmonics.The presented phenomena are in harmony with the deductions mentioned above.

A wide tube
Results for a tube with radius R 0 = 0.5 m are also presented.The calculations were made for frequency 100 kHz.
The cut-off frequency of the first symmetric mode for this tube is approximately equal to 911 Hz, thus the frequency of the described waves is considerably higher than the cut-off frequency.The Fermi distribution was chosen as the initial condition for the KZK equation.All other parameters of the computation were the same as in the previous case of the narrow tube.

Summary
This paper shows the derivation of the KZK equation in the parabolic approximation.The following validity condition of this equation for the waveguide has to be satisfied v ~m, v ^~m 2 (that is v v ^).The validity conditions of the KZK equation in the case of the description of quasi-linear waves propagation in the circular waveguide are further discussed by means of these conditions.It is shown that the condition kR 0 1 represents only the sufficient condition for validity of the KZK equation.The equation remains valid in the case of the small radius waveguide (and consequently with strong influence of the boundary layer) for frequencies below the cut-off frequency though kR 0 <1 under the condition that v v ^.The results of numerical solution of the KZK equation are presented for waveguides of both small and large radius which demonstrate these two cases.The numerical solution of the KZK equation was performed by means of a method based on the Bergen code [32], which was modified by the authors because it was necessary to incorporate the boundary condition, taking into account the boundary layer at the waveguide lateral walls.The new method for calculating the diffraction effects generated by the boundary layer effects is described in the text.When the plane wave with a frequency below the cut-off frequency propagates in a small radius tube, then its space evolution along the radius varies only a little and the results of the KZK and the Burgers equation are comparable.Solving the KZK equation also enables us to get the distribution of the perpendicular component of velocity.By means of the perpendicular component of the velocity we can see whether the validity limits of the KZK are fulfilled.In the waveguide of the large radius we can describe by means of the KZK equation the high frequency wave propagation for non-planar distribution of the velocity along the tube radius (it is no longer possible to use the generalized Burgers equation here).Again we can decide from the magnitude of the perpendicular component of the velocity whether we will use the KZK equation within the limits of its validity for the given parameters f, R 0 , v m and for the given velocity distribution along the radius near the source.

Fig. 8 :
Fig. 8: KZK equation solution of the longitudinal (V 3 ) and transversal (W 3 ) velocity component of the third harmonic