On the evolution equation for finite amplitude sound in relaxing media

Relaxation processes in a medium return thermodynamic parameters to equilibrium states. In the presence of acoustic waves, they can change the phase velocity and absorb acoustic energy. The problem of nonlinear propagation of finite amplitude sound through media with different relaxation processes had appeared within the scope of investigations more than half a century ago.^1–3 

The differential equation in particle velocity for the simplest case of a one-dimensional plane acoustic wave propagating in a medium with one basic relaxation mechanism was first derived in 1962 as follows:¹ 

where x and θ = t − x/c₀ are the space and time coordinates in a retarded frame, respectively, c₀ is the small-signal sound velocity, τ is the relaxation time, ε is the nonlinearity coefficient, the parameter $m=c∞2−c02/c02$⁠, c_∞ is the sound speed at infinite frequency in a medium with one relaxation process.

Equation (1) (or its modification with acoustic pressure p as the dependent variable³) was always considered to be unsolvable analytically except for some approximate or asymptotic cases;^1–3 therefore, at least one term was omitted for analysis of solutions, or a numerical approach was used.⁵ Solutions in the form of a transcendental Lambert W-function⁶ can also be considered as approximate. Thus, in the stationary case (∂v/∂x = 0), it was found that for relatively weak nonlinearity (⁠$D≡mc0/2εv≫1$⁠), Eq. (1) had a solution in the form of a shock with tanh-dependence on retarded time as follows:^1–3 

where v₀ was the acoustic velocity at the infinite value of θ, and another one boundary condition was dv/dθ → 0 at θ → ∞. It should be noted, as Rudenko and Soluyan admitted,² that the condition D ≫ 1, in fact, contradicted the assumption that the relaxation was small with respect to nonlinearity, which had been used in derivation of (1). In the opposite case (D ≪ 1), it is usually declared that (1) has an ambiguous solution, and some additional assumptions are considered and used.^1–3 

Since the 1960s, different new methods have been developed for solution of nonlinear equations of mathematical physics, and numerous books and papers have been devoted to that topic (a rather wide overview could be found in Refs. 7 and 8). It becomes possible to derive exact solutions to equations, which had been considered unsolvable analytically for decades, using not very complicated algebra.

Some attempts to find exact solutions to (1) are described in the current paper. It has been found that they were all unbounded; therefore, a modification of (1) is proposed so that the resultant equation could have solutions, interesting from the physical point of view.

As we had found by applying the classical Lie-group analysis,⁹ the equation under consideration admits only translation symmetries,

It is well known¹⁰ that application of conservation laws to differential equations can result in solutions other than invariant ones. Rather recent works by Ibragimov (Ref. 11 as an example) extended the range of equations, for which the conservation-law approach is applicable. Following that approach, we can formulate the conservation law for the first of operators, (3)

Thus, the next trivial exact solution of (1) can be easily found from

where $C̃1$ and $C̃2$ are the integration constants. Function (5) monotonously grows in time and diminishes in space.

By constructing the conservation law for the second of operators (3), we come to the next ordinary differential equation (Abel’s ODE of the second kind),

An absolutely identical equation was previously² derived from (1) in the stationary case (∂v/∂x = 0), and there is no contradiction because formulation of the conservation law for the operator X₂ from (3), in fact, means that the initial Eq. (1) is considered to be independent on the spatial variable x in this particular case. It was found in Ref. 2 that solution (2) could be obtained from the latter differential equation.

Applying non-classical symmetries¹² to (1), we have got another one set of transformation operators,

The system of characteristic equations corresponding to a linear combination of operators X₂ and X₃ looks like

The latter system yields two invariants,^9–11 the independent variable λ and the dependent Φ(λ), which are defined as follows:

Changing the variables in (1) with the help of (9) leads us to the next ordinary differential equation in function Φ,

By integrating this equation once, we arrive at the expression,

where K is the constant of integration. Expression (11) is the Abel equation of the second kind, which does not have explicit exact solutions in the general case.¹³ Nevertheless, as it is seen from its right-hand side, any function Φ(λ), if it ever satisfied (11), would tend to infinity at growing λ. The same conclusion is valid for approximate solutions of (11). Indeed, if nonlinear effects are much stronger than relaxation, then the second term on the left-hand side of (11) can be omitted, and the resultant ODE is solved analytically. For K = 0, we obtain the next exact solution of (11),

K₁ is another one integration constant. As, according to (9), acoustic velocity is proportional to the invariant variable Φ, it means that v is also unbounded.

However, we can try to solve Eq. (1) exactly in a much simpler way. Using the fact that it can be split into two equations, that is, linear and quadratic in v, we can search for a monochromatic harmonic solution,

where ω is the cyclic frequency of a wave, and w is its phase speed in the retarded frame.

When substituting expression (13) in (1), equalizing identical powers of harmonic functions and solving the obtained algebraic system of equations in amplitudes A and B, we can easily find two sets of exact solutions to (1),

and

where $j=−1$⁠.

Evidently, both expressions (14) and (15) are unbounded.

The fact that evolution Eq. (1) has nontrivial exact solutions, although they are not of interest from the physical point of view, suggests that maybe some revision of this equation could be desirable. For this purpose, we use the same Navier–Stokes equation as well as the continuity equation and the equation of state, as it was described in detail in Rudenko and Soluyan’s book,² 

to the same degree of approximation,

where ρ = ρ₀ + ρ′ is the mass density of a medium in the presence of propagating sound, ρ₀ is the density of a medium in equilibrium, κ is the thermal conductivity, c_v and c_p denote the specific heat at constant volume and pressure, respectively, μ_B and ζ are the bulk and shear coefficients of viscosity, respectively, p^′ is the acoustic pressure, and t′ is the integration variable. As it had been considered previously,^1–4 only the one-dimensional case is considered here for simplicity.

The same sequence of operations as in Ref. 2 is used, except for a change in the coordinate frame: ρ′ and p′ are excluded from the Navier–Stokes Eq. (16) using Eqs. (17) and (18) on condition (19). We do not neglect the thermal loss represented with the sound diffusivity δ, like it was also proposed previously.^2,14 Because it is known that usually $δ/(c02τ)≪1$ (look, for example, Ref. 4), terms of the second order in v are also kept on both sides of the final equation. Hence, we get

To get rid of the integral in the latter equation, it was differentiated by t and multiplied by τ; then, the result was summarized with (20). Thus, we get

Besides the coordinate frame, (21) differs from (1) in presence of the third and the fifth terms in the square brackets and of the last summand on the right-hand side. Elimination of the last term on the right-hand side was usually accounted for the same order of smallness of the parameter m as compared with any of the ratios in (19).^1,2 This restriction seems to be excessive.

Let us consider the outgoing solutions of (20): $v=fx−wt≡fχ$⁠. That leads to the next ODE in f(χ),

Coefficients M, N, P, and Q are determined as

After integration of the third-order ODE (22) we have

This equation does not pass the Kowalevskaya–Gambier test¹⁵ (the third step fails); it means that there is no general analytical solution of (24) although partial solutions can exist.⁸ 

For example, let us look for solutions of the next kind (outgoing weak shocks),

Substituting this expression in (24) and equating to zero coefficients before equal powers of exponential function, we obtain the next algebraic system,

Many different sub-cases are possible here. From the physical point of view, they correspond to different boundary conditions.

(a) If Δ = 0, then from the second of (26), α = −P, while from (29), K₁ = 0. From the sum of (27) and (28), it follows β = 4PQ/(2MQ + NP) so that the next quartic equation in w arises as

Even though the exact solutions of quartic equations are well known and can be written in the analytical form, the solution of the latter one is utterly cumbersome. For approximate calculation of roots, we neglect sound diffusivity (because $δ≪mc02τ≪c02τ$⁠) and then apply the perturbation method, considering m ≪ 1 (Ref. 4) as a small parameter. At the first step of perturbation, (30) is reduced to a cubic equation, which has three real roots: w₁ = c₀, w₂ = c₀(3 − 2ε), and w₃ = c₀(2ε − 1). The root w₂ is rejected because it has negative values for liquids (ε > 2). At the second step, we find parameters of solution (25) in the considered sub-case,

The last root w₃ is discarded because it leads to the like-sign α and β that contradicts the required solution in the form of outgoing shock. Approximate results (31) are compared with those calculated numerically from (30). We used material parameters for air at 20 °C: c₀ = 343 m/s, m = 1.25 × 10⁻⁴, δ = 1.86 × 10⁻⁵ m²/s, and τ = 4.73 × 10⁻⁴ s (relaxation time for N₂ molecules in air with 50% humidity).¹⁶ Dependence of normalized amplitude of acoustic speed f/mc₀ on χ = x − wt is shown in Fig. 1.

The solid curve corresponds to the numerical solution of (30), and the broken one corresponds to the approximate solution with parameters (31). It is seen that the height of a weak shock differs insignificantly, while the difference in the rise rate is more essential. Difference in w₁ values is less than 10⁻⁴%. Numerical calculation also showed that the fourth root of (30) is negative.

(b) If α = −Δ in the first of expressions (26), then for K₁ = 0, the variable w obeys the next algebraic equation,

Considering δ to be negligible, we arrive at the cubic equation in w, which has three real roots,

Hence, the next set of parameters of solution (25) can be found in this sub-case for the only root, which remains positive for all values of ε,

(c) The integration constant K₁ in (24) is not obligatorily equal to zero. Suppose that, for example,

where r is a positive real number. In this case, if α = − Δ, the amplitude of a shock obeys

and if additionally |4QK₁| ≪ P², then we have

Finally, we can conclude that if α and β have opposite signs and w is positive, then any solution given by (25) with parameters obeying (31) or (33)–(37) corresponds to an outgoing weak shock. However, it should be mentioned that (11) allows solutions in form of a hyperbolic tangent too,

When repeating similar algebra, we find that $α=±K1−2c02ε−12/m1/2$⁠. If, for example, $K1=6c02ε−12/m$⁠, then we obtain rather simple analytic expressions,

for

A wide range of integration constantsK₁ is also possible when searching for solutions of (24) in the form of a hyperbolic tangent (38) that allows different boundary conditions.

Examination of consequences of the evolution equation for finite amplitude acoustic waves in a relaxing medium suggested correction of that equation. It was shown that taking into account the thermal and viscous loss and the usually omitted nonlinear terms of the second order leads to a new third order partial differential equation, which has some nontrivial analytical solutions, at least solutions in the form of weak shocks.

The data that support the findings of this study are available within the article.