Ultrasonic Waves in Bubbly Liquids: An Analytic Approach

: We consider the problem of the propagation of high-intensity acoustic waves in a bubble layer consisting of spherical bubbles of identical size with a uniform distribution. The mathematical model is a coupled system of partial differential equations for the acoustic pressure and the instantaneous radius of the bubbles consisting of the wave equation coupled with the Rayleigh–Plesset equation. We perform an analytic analysis based on the study of Lie symmetries for this system of equations, concentrating our attention on the traveling wave case. We then consider mappings of the resulting reductions onto equations deﬁning elliptic functions, and special cases thereof, for example, solvable in terms of hyperbolic functions. In this way, we construct exact solutions of the system of partial differential equations under consideration. We believe this to be the ﬁrst analytic study of this particular mathematical model.


Introduction: Mathematical Model
We consider an analytic approach to the study of high-intensity nonlinear ultrasonic waves in bubbly liquids, a problem which has many applications in, for example, medicine, industry, and particle filtration. The known mathematical models to describe this problem consist of a coupled system of partial differential equations [1]. The first is the linear nondissipative wave equation which describes acoustic pressure changes in the bubble layer (see, for example, [1][2][3][4][5][6][7]). The second equation, which governs the bubble radius changes, is the Rayleigh-Plesset equation [8][9][10], which has been generalized in various ways and is sometimes written in terms of the bubble volume variations (see, for example, [3,5,7,11]). As can thus be seen, this problem is one that has attracted much interest over the years, due to its important practical interest, with many studies having been undertaken from a numerical point of view.
In this paper, of the various models that have been proposed, we follow most closely that in [1]. We assume a liquid layer of spherical bubbles of identical size distributed uniformly, and consider the interaction between the acoustic waves and the bubble oscillations. The sound speed inside the bubble layer is denoted by c L . We assume that the density in the layer can be approximated by the density outside the layer ρ 0 because of small differences between the density of water at the equilibrium state and in the bubbly liquid.
The acoustic pressure P = P(x, t) in the bubble layer obeys the linear wave equation [1][2][3][4][5][6][7] where β = β(x, t) is the local fraction of the volume occupied by the gas. If we now assume a constant number N of air bubbles per unit volume, the function β is given by R = R(x, t) being the instantaneous radius of the bubbles. The wave equation can then be written as We now consider the nonlinear interaction of ultrasonic waves with the bubbles in the layer. The mathematical model describing this interaction is given by the wave Equation (3) coupled with the Rayleigh-Plesset equation written in terms of the instantaneous radius of the bubbles [1] RR tt + 3 2 where p g = 2 σ R 0 + p s − p v . In the above equation, p v is the gas pressure inside a bubble, p s is the ambient static pressure, R 0 is the equilibrium bubble radius, γ is the polytropic exponent of the gas, and σ is the coefficient of surface tension. We note that we do not include a viscosity term in Equation (4). We also note that the physically possible values of the polytropic exponent γ are 1 ≤ γ ≤ 5/3, where the case γ = 1 corresponds to the isothermal case and the case γ = 5/3 to the adiabatic case.
The model is therefore given by the system of coupled partial differential Equations (3) and (4) for the acoustic pressure P(x, t) and the instantaneous radius of the bubbles R(x, t), both of which are considered to be functions of space and time coordinates x and t. The aim of this paper is to study this model from an analytic point of view. Whilst various aspects of this model have been explored in the literature, we believe the work presented here to constitute the first analytic results for this system of partial differential equations. Similar approaches have been used to obtain analytical or semi-analytical solutions of physically interesting problems in [12,13]. In Section 2, we perform a Lie symmetry analysis of this system of equations and obtain the possible similarity reductions. In Section 3, we analyze in detail the traveling wave reduction, for which we obtain solutions by mapping onto equations defining elliptic functions (see [14,15]), and special cases thereof, for example, solvable in terms of hyperbolic (or trigonometric) functions; lists of the transformed equations are given in the Appendices A and B. In this way, we obtain exact solutions of the system of partial differential Equations (3) and (4). In the final section, we discuss the results obtained.

Similarity Reductions
We now consider the derivation of similarity reductions associated with classical Lie point symmetries [16][17][18] for the system of Equations (3) and (4). In order to do so, we require the invariance of this system under the one-parameter Lie group of infinitesimal transformations in (x, t, P, R) given by where is the group parameter. The symmetry generator associated with the above group of point transformations can be written as The invariance condition leads to an overdetermined system of linear differential equations (the determining equations) for the infinitesimals ξ, τ, φ 1 and φ 2 . Once the in-finitesimals have been obtained, the similarity variables are found by solving the associated characteristic equations dx ξ(x, t, P, R) The infinitesimals ξ, τ, φ 1 , and φ 2 associated with the classical Lie symmetries of the system of Equations (3) and (4) depend on the parameters appearing in the equations. In the general case, in which there is no restriction on the physical parameters, we find that the infinitesimals are where c 1 and c 2 are arbitrary constants; this implies that, when solving the characteristic Equations (10), the only possible similarity reduction is the traveling wave reduction. We will analyze this reduction in full detail in Section 3.
In the particular case in which p v = p s and σ = 0, which implies p g = 0, we obtain as infinitesimals which provide, in addition to the traveling wave reduction, a scaling reduction. This scaling reduction corresponds to the choice c 3 = 0, in which case we can set c 1 = c 2 = 0 and c 3 = 1 without lost of generality. The similarity variables are and the associated system of ordinary differential equations is where primes denote derivatives with respect to z. Solving Equation (21) for p and substituting in Equation (20) yields a fourth order ordinary differential equation for r. This could be solved numerically, for appropriate initial conditions, or at least for a half-set of initial conditions with optimal values for others to be determined via numerical experiments. Alternatively, the change of independent variable z = e y allows us to write this system as ρ 0 πN(42r 3 + 6rr 2 + 3r 2 r + 15r 2 r ) = 0, (22) p + ρ 0 (12r 2 + 3rr + rr + 3 2 where now prime denotes derivatives with respect to y. One possibility for obtaining semianalytical solutions of this system would be to consider when y ∼ 0 or y → ∞. For both of these cases, solving (23) for p and substituting in (22) yields an autonomous fourth order ordinary differential equation for r. Being autonomous, we can reduce its order using the change of variables r (y) = f (r); this then yields a third order ordinary differential equation for f (r). We will return to a study of the system (22), (23) in a later paper.

Travelling Wave Reduction: Exact Solutions
In this section, we consider the traveling wave reduction for the system of coupled partial differential Equations (3) and (4). Under the similarity reduction we obtain the system of ordinary differential equations where primes denote derivatives with respect to z and the constants C, a 0 , a 1 and a 2 are defined in terms of the physical parameters via We note here that, in the special case p s = p v and σ = 0, which implies p g = 0, we have a 0 = a 1 = a 2 = 0. We also note that, in (27), the wave speed c is assumed to be different from the sound speed inside the bubble layer c L .
The first equation can be integrated twice to give where A and B are the two arbitrary constants of integration. Inserting the above expression for p into the second equation, we obtain the second order ordinary differential equation for r rr + 3 2 The study of this equation splits into the two cases γ = 1 (the isothermal case) and γ = 1. Below, for each of these two cases, and with A = 0, we integrate this equation to obtain Equations (34) and (47), respectively. Let us note that, alternatively, we could in this autonomous case A = 0 use the change of variable r (z) = F(r) in order to obtain the Bernoulli equation The transformation G = F 2 yields a linear equation for G. Solving for G and then inverting the changes of variable G = F 2 = (r (z)) 2 leads again to Equations (34) and (47).
We now consider first the case γ = 1.

The Non-Isothermal Case γ = 1
We consider the autonomous case A = 0 of Equation (32). Multiplying by 2r 2 r and integrating, we obtain where D is the constant of integration. Here, we have assumed that γ = 1, i.e., we have excluded the isothermal case. We will come back to the isothermal case in the next subsection.
In the case γ = 1, we then have that the variable r satisfies the first order Equation (34). We now rescale r via to obtain where the constants B 1 , B 2 , B 3 and B 4 are defined in terms of the physical parameters a 0 , a 1 , a 2 and C and the constants of integration B and D as (40) Our objective now is to find solutions for the system of Equations (3) and (4) by searching for solutions of its traveling wave reduction, which in the autonomous case A = 0 is given by Equation (36). Equation (36), or equivalently (34), allows us to obtain such solutions via a quadrature. Here, instead, we search for transformations of (36) to ordinary differential equations with known general solutions. In particular, we look for transformations that map Equation (36) onto a first order nonlinear ordinary differential equation of the second degree of the form in new dependent and independent variables V = V(τ), where P (V) is a polynomial in V of degree less than or equal to four and primes denote derivatives with respect to the new independent variable τ (for the general case, i.e., in order to give rise to elliptic functions, this polynomial needs to be of degree four or three; when this polynomial is of degree two or one, the equation is readily solvable in terms of elementary functions). In order to do so, we consider, as in [19][20][21], the combined Sundman power-type transformation where δ and ε are real numbers and ε = 0 (we note that the Sundman transformation was originally proposed in [22]; see also [23]). Applying the transformation (42) to Equation (36) yields where µ = 2(1 + εδ). We then require that the right-hand side of Equation (43) is a polynomial in V of degree at most four. We may then express solutions using Jacobi elliptic functions, Weierstrass elliptic functions, hyperbolic (or trigonometric) functions, or even as rational functions. The analysis splits into two different cases depending on whether ε is positive or negative, since the sign of ε determines whether the powers of V in (43) are decreasing or increasing. However, without loss of generality, the case ε < 0 may be discarded, as we now explain. We set m = ε + µ to obtain If we now set V = 1/U, this equation then becomes where n = 4 − m. Equations (44) and (45) are of the same form but with opposite signs of ε. We see that Equations (43) of the required form (right-hand side polynomial in V of degree at most four) for ε < 0 can be obtained from those having ε > 0 via the transformation V = 1/U. Moreover, given the relation V = 1/U between the solutions of these equations, and since δ is left invariant, we see that the only difference in the solutions of (36) obtained using the transformation (42) is in the sign of z. Since Equations (34) and (36) are invariant under z → −z, this means that the solutions obtained in the cases ε > 0 and ε < 0 are trivially related. For reasons of clarity, we leave the complete classification of Equations (43) having a right-hand side polynomial in V of a degree of at most four, in the case ε > 0, to Appendix A. We note that, in this classification process, we may set B 2 = 0 or B 4 = 0 by an appropriate choice of the constants of integration B and D. However, setting B 1 = 0 implies that the coefficient of surface tension σ is zero, whereas setting B 3 = 0 means that

The Isothermal Case γ = 1
We now turn our attention back to the isothermal case, that is, the case where the polytropic exponent γ = 1. Considering once again the autonomous case A = 0 of Equation (32), we see that the variable r satisfies As before we multiply by 2r 2 r and integrate, this time obtaining where E is the constant of integration. This equation then allows us to obtain solutions of the traveling wave reduction of the system (3) and (4) with γ = 1, in the autonomous case A = 0, via a quadrature. Here, instead, we search for transformations to ordinary differential equations with known general solutions. Setting a 2 = 0, in order to remove the logarithmic term, we obtain the equation This can be compared to Equation (34) in the special case a 2 = 0. As previously, we rescale r via to obtain where the constants B 1 , B 2 , and B 4 are defined as that is, corresponding to Equation (36) with B 3 = 0 (and D replaced by E). Using once again the combined Sundman power-type transformation (42), we obtain where µ = 2(1 + εδ), i.e., corresponding to Equation (43) with B 3 = 0. Again as before, we need only consider the case ε > 0. The classification of equations of the form (54), with a right-hand side a polynomial in V of a degree of at most four, in the case ε > 0 is described in Appendix B. In the same way as in the previous subsection, our aim is to use such equations in order to obtain solutions of the system of Equations (3) and (4), now in the special case where γ = 1 but also a 2 = 0, i.e., p g = 2 σ R 0 + p s − p v = 0. Comparing with the previous subsection, we see that the equations resulting from this classification are in fact valid for a 2 = 0 and any 1 ≤ γ ≤ 5/3. We recall that the case γ = 1 and a 2 = 0 is governed by Equation (47).

Examples
We now give examples of the exact solutions that it is possible to obtain from the equations given in Appendices A and B. We present examples of exact solutions given in terms of elliptic functions, as well as in terms of hyperbolic and trigonometric functions. For these last cases, when the right-hand side of (43) is a cubic or quartic polynomial in V, we may need to impose conditions on the coefficients B 1 , B 2 , B 3 and B 4 of the equations in Appendix A in order for this polynomial to have a double root. When imposing any conditions on these four coefficients, we note that we are free to choose B 2 and B 4 , as these depend on the constants of integration B and D/E. However, B 1 and B 3 depend on physical parameters, although their dependence on the wave speed c may allow greater freedom in imposing conditions involving these two constants. Again, we remark that setting B 1 = 0 implies σ = 0, and setting B 3 = 0 implies p g = 2 σ R 0 + p s − p v = 0.

Example 1
We first consider Equation (A7) given in Appendix A, for which ε = 1/2 and δ = 3/2, and where, since we assume B 3 = 0, we have γ = 5/3, i.e., our choice of example corresponds to the adiabatic case. In the general case, when the quartic polynomial on the right-hand side of Equation (55) has four distinct roots, the general solution can be expressed using Jacobi elliptic functions. Alternatively, we may express the general solution in terms of Weierstrass elliptic functions, as where Here, we focus on the case where the polynomial in V on the right-hand side of (55) has a double root, i.e., we assume that the equation is of the form Comparing coefficients, we see that we must have c = −2a − b, so our equation in fact reads and The first two of these equations yield and thus we obtain a correspondence between (B 1 , B 3 ) and (a, b). The third of the above equations can then always be satisfied through a choice of the constant of integration D on which B 4 depends. We now assume B 1 < 0, B 3 < 0 (for γ = 5/3, they have the same sign) and B 2 1 + 12B 3 > 0. It then follows that both a 2 and (a + b) 2 as obtained above are real and positive, and so we may solve these relations for real a and b. Equation (59) has a general solution where τ 0 is the constant of integration. We note that and so choosing upper signs in (63) and (64) yields a solution defined in terms of hyperbolic functions, whereas choosing lower signs will give a solution defined in terms of trigonometric functions. We note that, for a choice of upper signs, we have and so, in order to avoid singularities, we solve (63) and (64) for a and (a + b) of opposite sign. We thus obtain (here ± corresponds to the sign of a), where Solutions having singularities are given by (again, ± corresponds to the sign of a).
Choosing lower signs in (63) and (64), we have and so singularities cannot be avoided. We obtain for any a, (a + b) defined by (63) and (64). We note that, for the upper choice of sign in (68), V(τ) as given therein is always positive. We may then write a solution u(z) parametrically as We then find r(z) and p(z) from (35) and (31) respectively (with A = 0), and thus the corresponding traveling wave solution (24) of the system (3) and (4).

Example 2
We now consider Equation (A25) given in Appendix A, for which ε = 1/4 and δ = −1/2, and where we assume that B 3 = 0 and therefore γ = 5/3. Assuming B 1 < 0 and B 3 < 0 (for γ = 5/3 they have the same sign), we obtain the general solution of this equation as where τ 0 is the constant of integration. We note that V(τ) is always positive, and we may therefore write a solution u(z) parametrically as We then find r(z) and p(z) from (35) and (31) respectively (with A = 0), and thus the corresponding traveling wave solution (24) of the system (3) and (4).

Example 3
We now consider Equation (A23), for which ε = 1/3 and δ = −1/2. This equation belongs to both Appendices A and B, i.e., it corresponds to any value of the polytropic exponent γ in the range 1 ≤ γ ≤ 5/3. Equation (77) has the same form as Equation (74), but with B 1 and B 3 replaced by B 2 and B 4 , respectively. We choose the constant of integration in B 4 (i.e., D or E) so that B 4 < 0, and thus obtain the general solution of (77) as where τ 0 is the constant of integration. Since V(τ) is always positive, we may then write a solution u(z) parametrically as We then find r(z) and p(z) from (35) and (31), respectively (with A = 0), and thus the corresponding traveling wave solution (24) of the system (3) and (4). This solution is for B 1 = 0 and B 3 = 0, which corresponds to σ = 0, i.e., zero coefficient of surface tension, and p s = p v , i.e., gas pressure inside the bubbles equal to the ambient static pressure.

Example 4
We now consider Equation (A15), for which ε = 1/3 and δ = 1. This equation belongs to both Appendices A and B, i.e., it corresponds to any value of the polytropic exponent γ in the range 1 ≤ γ ≤ 5/3. The general solution of Equation (80) is given by where Here, we focus our attention on the case where the polynomial in V on the right-hand side of (80) has a nonzero double root and so we compare with which implies and therefore the coefficients B 2 and B 4 are subject to the constraint B 2 2 = 4B 4 , a relation that can be satisfied by an appropriate choice of the constant of integration in B 4 (i.e., D/E) in terms of that in B 2 (i.e., B). We assume B 2 < 0 (through an appropriate choice of the constant of integration B) and thus obtain the general solution of (80) as where τ 0 is the constant of integration. This solution V(τ) is always positive, and we may write a solution u(z) parametrically as We then find r(z) and p(z) from (35) and (31), respectively (with A = 0), and thus the corresponding traveling wave solution (24) of the system (3) and (4). This solution is for B 1 = 0 and B 3 = 0, which corresponds to σ = 0, i.e., zero coefficient of surface tension, and p s = p v , i.e., gas pressure inside the bubbles equal to the ambient static pressure.

Example 5
We finally consider Equation (A30), for which ε = 1/4 and δ = −5/2. This equation also belongs to both Appendices A and B, corresponding to any value of the polytropic exponent γ in the range 1 ≤ γ ≤ 5/3. The general solution of (88) is where τ 0 is the constant of integration. Assuming now that B 1 < 0, V(τ) is always positive, and we can thus write a solution for u(z) parametrically as We then find r(z) and p(z) from (35) and (31) respectively (with A = 0), and thus the corresponding traveling wave solution (24) of the system (3) and (4). This solution is for B 3 = 0, which corresponds to p g = 2 σ R 0 + p s − p v = 0.

Discussion
We have undertaken a study of a system of partial differential equations describing the propagation of high-intensity acoustic waves in a bubble layer. We believe our work here, which has allowed the derivation of several exact solutions, to be the first analytic study of this model. In later papers, we will further explore some of the examples considered here, as well as use other examples of transformed equations to derive further exact solutions. We will also return in later papers to the consideration of the system (22) and (23), which we obtained as a second similarity reduction of the system (3) and (4), as well as of the non-autonomous case A = 0 of Equation (32). Within this context, but also more generally, it would be of interest to explore transformations to other classes of solvable/integrable differential equations.
Author Contributions: Conceptualization, P.R.G. and A.P.; methodology, P.R.G. and A.P.; investigation, P.R.G. and A.P.; writing-original draft preparation, P.R.G. and A.P.; writing-review and editing, P.R.G. and A.P. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
We thank the referees for their helpful remarks and kind suggestions, in particular with regard to semi-analytical procedures.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A. Classification: γ = 1 Here, we give the results of our classification of equations of the form (43) in the case ε > 0. We first recall that the parameter δ in the transformation (42) is given by δ = (µ − 2)/(2ε). Since the possible values of the polytropic exponent γ are 1 < γ ≤ 5/3 (the case γ = 1 is considered separately in Appendix B), we see that the terms on the right-hand side of Equation (43) are written in order of (strictly) decreasing powers of V. We therefore take our equation in the form (43), and consider the different possible cases where the right-hand side is a polynomial in V of degree 4, 3, 2 and 1. In the resulting four lists of equations presented below, the first coefficient B i , i = 2, 1, 4 to appear on the right-hand side is assumed to be nonzero (by the first coefficient, we mean that which is furthest to the left). When the first coefficient on the right-hand side is B 3 , which occurs for the last equation in each list, this coefficient may be zero or nonzero.
First of all, we consider the case in which the right-hand side of Equation (A1) is a quartic polynomial in V, which corresponds to the choice µ = 4 − ε, which then implies δ = (2 − ε)/(2ε). We obtain the following equations: For each of the above (except the last with B 3 = 0), the value of ε is fixed. For three of the above equations, if B 3 = 0, then γ is specified as γ = 5/3. For B 3 = 0, these three equations are valid for any polytropic exponent 1 < γ ≤ 5/3, as is the case for all of the other equations.
We now consider the case in which the right-hand side of Equation (A1) is a cubic polynomial in V, which corresponds to the choice µ = 3 − ε, which then implies δ = (1 − ε)/(2ε). We obtain the following equations: The same considerations as for the quartic polynomial case apply here. For each of the above (except the last with B 3 = 0), the value of ε is fixed. For Equation (A18), if B 3 = 0, then γ is specified as γ = 5/3. For B 3 = 0, this equation is valid for any polytropic exponent 1 < γ ≤ 5/3, as is the case for all of the other equations.
In the case in which the right-hand side of Equation (A1) is a quadratic polynomial in V, which corresponds to the choice µ = 2 − ε, and which then implies δ = −1/2, we obtain the following set of equations for V: As in the quartic and cubic cases, for each of the above (except the last with B 3 = 0), the value of ε is fixed. For Equation (A25), if B 3 = 0, then γ is specified as γ = 5/3. For B 3 = 0, this equation is valid for any polytropic exponent 1 < γ ≤ 5/3, as is the case for all of the other equations.
For the case where the right-hand side of Equation (A1) is a linear polynomial in V, which corresponds to the choice µ = 1 − ε, and which then implies δ = −(ε + 1)/(2ε), we obtain: As for the previous three classes of polynomial considered above, for each of the above equations (except the last with B 3 = 0), the value of ε is fixed. All of the above are valid for any polytropic exponent 1 < γ ≤ 5/3.

Appendix B. Classification: γ = 1
In this second Appendix, we briefly describe the results of our classification of equations of the form (54) in the case ε > 0. Again, we recall that the parameter δ in the transformation (42) is given by δ = (µ − 2)/(2ε). The polytropic exponent is γ = 1. The terms on the right-hand side of Equation (54) are written in order of (strictly) descreasing powers of V. We therefore take our equation in the form (54), and consider the different possible cases where the right-hand side is a polynomial in V of degrees 4, 3, 2, and 1. A comparison of Equation (A33) with Equation (A1) leads to the result that the sought-after equations in this isothermal case γ = 1 are precisely those given in Appendix A in the special case B 3 = 0. That is, our list of equations consists of all equations in Appendix A, where B 3 does not appear as a coefficient, as well as those resulting from setting B 3 = 0 where it does appear as a coefficient, always with γ = 1. However, since the choice B 3 = 0 can always be made in Appendix A, the resulting equations are in fact valid for any 1 ≤ γ ≤ 5/3.