Linear Iterative Procedure to Solve a Rayleigh–Plesset Equation

: A nonlinear Rayleigh–Plesset equation for describing the behavior of a gas bubble in an acoustic ﬁeld written in terms of bubble-volume variation is solved through a linear iterative procedure. The model is validated, and its accuracy and fast convergence are shown through the analysis of several examples of different physical meanings. The simplicity and usefulness of the presented method here in relation to the direct resolution of the whole nonlinear system, which is also discussed, make the method very attractive to solving a problem. This iterative method allows us to solve only linear systems instead of the nonlinear differential problem. Moreover, the implementation of the iterative algorithm includes a tolerance-dependent stopping criterion that is also tested.


Introduction
Rayleigh-Plesset (RP) equations model the behavior of gas bubbles impinged by an acoustic wave under the consideration of some physical approximations [1,2]. Two frameworks exist in this context. On the one hand, the radius model uses time-dependent bubble radius R(t) in the derivation of an ordinary differential equation. On the other hand, the volume model uses bubble volume variation v(t) as a time-dependent variable for deriving an ordinary differential equation (RPv). In this work, we consider the latest. Such nonlinear differential model is of second-order. Two initial conditions on primary variable v and its first time derivative v are required at the onset of the physical phenomenon (t = 0). This model equation was derived in 1973 by Zabolotskaya et al. [3]. It was used in conjunction with the wave equation to analyze effects caused by a population of gas bubbles in a liquid, which parameter of nonlinearity increases by several orders of magnitude, on an ultrasonic field through the development of second-order perturbation methods [3,4] and numerical models [5][6][7][8].
The iterative procedure applied in this work to solve the RPv equation follows the iterative method developed for solving problems defined by nonlinear second-order ordinary differential equations [9]. This method relies on an iterative procedure that considers the differential equation at each iteration by assuming its nonlinear terms as known from the previous iteration, which leads to the expression of a modified reformulated linear problem at each iteration, so that only linear differential equations must be solved throughout the entire process (no nonlinear ordinary differential equations have to be explicitly solved). This latter point is an advantage regarding other techniques that need to solve nonlinear systems at each iteration step. When integrals at each iteration step can be analytically evaluated or by means of symbolic computations, it allows for the solution to be expressed in polynomial form. When these integrals cannot be analytically or symbolically solved (such as in this work, see Section 2), an alternative method consists of numerically calculating the integrals.
In this paper, we develop an iterative procedure to solve the initial-value Rayleigh-Plesset problem defined by the RPv equation to track in time the oscillations of a single bubble in an ultrasonic field. Section 2 describes the application of the technique. To this end, the system obtained at each iteration is numerically solved. Analysis of several physical examples and a comparison to the reference data are presented in Section 3, thus showing its fast convergence and accuracy. The simplicity and usefulness of the method are also discussed through these examples. A stopping criterion is also developed in the implemented algorithm. Lastly, Section 4 concludes this work.

Materials and Methods
We consider a tiny air bubble placed in liquid and impinged by an ultrasonic field. To describe the response of the bubble to the acoustic field, we work within the Rayleigh-Plesset framework. In this context, some approximations are applied to model the behavior of bubble oscillations, and several physical assumptions are made: • bubble content is gas without vapor; • the bubble is spherical (spherically symmetric pulsation); • the gas contained in the bubble has adiabatic behavior; • surface tension at the gas-liquid boundary is neglected; • a quadratic Taylor expansion of the adiabatic gas law limits bubble oscillations to moderate values; • the bubble does not itself radiate sound.
Moreover, in this work, we neglect buoyancy, Bjerknes, and viscous drag forces. The evolution with time t of the volume variation of the bubble v(t) = V b (t) − V b0 , a smooth real-value function of t, where V b and V b0 = 4πR 3 b0 /3 are the current and initial volumes of the bubble, respectively, and R b0 is the initial radius, is modeled by the Rayleigh-Plesset equation (RPv) [3,4], in which v and v are the first-and second-order derivatives of v(t), is the resonance frequency of the bubble, where ρ l is the density of the liquid at the equilibrium state, γ is the specific heat ratio of the gas, p b = ρ b c 2 b /γ is its atmospheric pressure, ρ b and c b are the density and sound speed at the equilibrium state of the gas, and δ = 4ν l /ω b R 2 b0 is the viscous damping coefficient, for which ν l is the kinematic viscosity of the liquid. a = (γ + 1)ω 2 b /2V b0 is the second-order nonlinear parameter due to the adiabatic law used to model the gas pressure inside the bubble. b = 1/6V b0 is the nonlinear coefficient associated to the bubble dynamics. The parameter η = 4πR b0 /ρ l is a constant that affects the source term p(t), which is the time-dependent acoustic pressure perturbation that impinges the bubble (a known source function). Several source functions are considered in Section 3. t = T is the upper limit of the time interval. The derivation of this ordinary differential equation, Equation (1), can be found in [3,4].
The differential system is closed by imposing initial conditions expressing the rest of the bubble at the onset of the study (nonperturbation of the bubble at the start of the experiment), The iterative procedure relies on the resolution of several successive systems DS n , which include the linear terms in the ordinary differential equation (linear operator L), the terms in the ordinary differential equation that do not depend on the dependent variable v (function g), and the initial conditions of the differential system (auxiliary operator A), through an n−iterative process that starts at n = 0 without considering nonlinear terms in the ordinary differential equation (nonlinear operator N), giving the solution v 0 (t), and follows for n = 1, . . . , N by taking into account the nonlinear operator N evaluated at the preceding iteration, giving the solution v n (t): DS n (n = 1, . . . , N); Its application to Equations (1) and (2) requires the definition of linear operator L (linear terms of the differential equation), nonlinear operator N (nonlinear terms), known function g (independent terms on v), and auxiliary operator A (initial conditions), respectively: In this case, successive iterations are, for DS 0 , and for DS n , Each differential system obtained at each iteration n, Equations (9) and (10), is defined by a linear (in terms of v 0 and v n , respectively) inhomogeneous second-order ordinary differential equation of time-dependent variable v 0 and v n , respectively, and Cauchy conditions. Finding the analytical solution of Equation (10) is not straightforward. Symbolic calculus was used to obtain this solution at each iteration of the progression, but it failed in giving v n from n = 2 due to the complexity of the second-hand side of Equation (10).
Thus, we chose to apply a numerical scheme at each iteration n, namely, the secondorder explicit vector Euler approximation method [10], which is more efficient than solving the whole problem directly, Equations (1) and (2), numerically. A fourth-order vector Runge-Kutta algorithm was also developed, but the obtained results in the following section were not affected by the choice of numerical method. To this purpose, we reduce DS 0 and DS n to a nonlinear first-order system of two coupled ordinary differential equations by setting the new variable v (t) = V(t), which yields These differential systems are written in the following vector form, after defining The application of the above-mentioned numerical algorithm to Equations (13) (13) and (14) is sought at all points t i of that set (subindex i indicates the approximate value of the corresponding variable at point t i ). This process leads to the following formulation: The derivative V n−1,i in Equation (16), V 0,i when n = 1, is evaluated at each point t i through a first-order finite-difference formula [10].
Since no exact solution of our differential problem is known, we evaluate the error produced during the iterative process through the discrete L 2 -norm of the difference between two consecutive approximate solutions, V n (i, 1) and V n−1 (i, 1) on the one hand, i.e., v n,i and v n−1,i , and between V n (i, 2) and V n−1 (i, 2) on the other hand, i.e., v n,i and v n−1,i , at each point i of the discretization set t i for i = 1, . . . , I, and for all values n = 1, . . . , N, which is the E L 2 V n -vector of components: The objective of this paper is to propose a very useful tool for solving the problem described above. In this context, a stopping criterion is set from the value of the tolerance, , as follows. Iterations stop as soon as both inequalities stand: for which tolerance could be defined as, for instance, = v b0 × 10 −10 . An example of using the stopping criterion is given in Section 3. The ad hoc implementation of the scheme is performed in a MATLAB R2017a environment. The schematic workflow for the successive steps of the entire algorithm is represented in Figure 1. CPU times required to obtain the solution of the bubble problems from the resolution process, indicated by t cpu in the following paragraphs, are obtained by running the code on a 128 GB RAM computer with a 64 bit Windows 10 operative system working with an Intel Core™ i7-6800 K CPU @ 3.40 GHz.

Results and Discussion
In this section, the method developed in Section 2 is applied to solve the following physical problem: in water (ρ l = 1000 kg/m 3 , ν l = 1.43 × 10 −6 m 2 /s), an air bubble of initial radius R b0  Figure 2d,e, respectively, for v n and v n . The convergence of the method is fast (Figure 2d,e). Figure 3 displays the comparison of the obtained results here to the obtained data by a simulation carried out with the Snow-Bl code considering one single bubble with the same parameters [5]. This code relies on the approximation of the solution of Equation (1), coupled to the linear wave equation, by finite-difference approximations. It can track both pressure and bubble oscillation fields in time. The concordance of both bubble-volume variation results validates the method developed in this paper. The CPU time needed by the Snow-Bl code is 4893 s, which is 12 times the time required here.   Figure 4 shows the corresponding results for the infinitesimal source amplitude p 0 = 1 × 10 −3 Pa, Case 1 , which implies a linear behavior of the bubble response since the nonlinear parameters a and b are multiplied by very small values. Convergence is thus immediate, after one iteration. Figure 5 shows the corresponding results for the same finite source amplitude as that in Case 1, p 0 = 1.25 kPa, but after setting nonlinear parameters a and b to zero, Case 1 . This imposition implies a linear response of the bubble even at this high source amplitude since ordinary differential Equation (1) is linear. The solution is reached at the first calculation for v 0 and v 0 , and no iterations are needed. The comparison of the solution, especially v, for Cases 1 and 1 reveals that both curves are identical (proportionally to their amplitude) in Figures 4 and 5, i.e., the result is similar in both cases. However, in Figure 2 for Case 1, the nonlinearity due to finite amplitude and nonlinear parameters a and b affects the system and modifies the solution iteration by iteration, which leads to a different signal, showing distortion and a slight nonsymmetrical response between positive and negative volume variations around axis v = 0.  To illustrate the possibilities and usefulness of the method implemented in this paper, we consider two other forcing terms in Equation (1)   For Case 1, we test the stopping criterion given by (18) with = v b0 × 10 −1 , = v b0 × 10 −10 , and = v b0 × 10 −20 . Results are shown in Figure 8. N = 7 with t c pu = 105 s, N = 26 with t cpu = 1424 s, and N = 63 with t cpu = 8450 s are respectively required for satisfying the tolerance. Figure 8c shows an extremely exigent case since it takes many iterations to slightly decrease error and improve results v n and v n . In the three cases presented in this section, even though the parameters of the problem set two very high nonlinear coefficients, a = 6.95 × 10 +28 Hz 2 /m 3 and b = 4.37 × 10 +14 /m 3 , the convergence of the presented method is fast, indicating the good behavior of the method. Regarding the benefits obtained here to solve the RPv in relation to other potential or existing resolution methods, the following arguments are given. (i) First, the direct treatment of the complete nonlinear differential system by reducing its order to 1 would end up requiring a numerical resolution (vector Newton-Raphson or similar algorithm) of a nonlinear coupled system of two discretized equations (for which the nonlinear operator N would depend on the unknown variable N(v n ) instead of the known solution at the previous iteration N(v n−1 ) in Equation (10) of the current procedure) at each step of the numerical process (vector Runge-Kutta or vector explicit Euler method), of which the time cost is elevated, and accuracy depends on the choice of the required vector estimation to start the system resolution. A clear advantage of the presented model is that it does not need such intermediate solution to nonlinear systems since the nonlinear terms used are known from the previous iteration. (ii) On the other hand, the direct treatment of the complete nonlinear differential system of second order, Equations (1) and (2), is much likely to be more delicate and time-consuming. The direct substitution of derivatives by finite-difference formulas in Equations (1) and (2) would lead to an implicit scheme requiring the resolution of a nonlinear system of equations at each time step of the process, which is usually highly time-consuming. The treatment of differential system Equations (1) and (2) by finite-volume time approximations would end up to a quite similar resolution of nonlinear systems of equations at the time steps. The use of finite-element approximations through the application of a variational formulation to obtain a weighted-integral form equivalent to Equations (1) and (2) would also require the solution of a nonlinear system of equations at each step of the time discretization process (finite differences). Thus, the simplicity of the procedure presented in this paper makes it very useful and attractive.

Conclusions
A nonlinear Rayleigh-Plesset equation written in the bubble variation volume framework was solved through a linear iterative procedure. The presented results showed fast convergence and accuracy of the technique developed in this paper. The simplicity and usefulness of the method presented here in relation to the direct resolution of the whole nonlinear system was also discussed. This iterative method allowed us to solve only linear systems instead of a nonlinear differential problem. Moreover, the implementation of the iterative algorithm includes a tolerance-dependent stopping criterion that was also tested.