Transient effects in the translation of bubbles insonated with acoustic pulses of finite duration

# Transient effects in the translation of bubbles insonated with acoustic pulses of finite duration

Patricia Vega-Martínez1
Javier Rodríguez-Rodríguez1 Email address for correspondence: javier.rodriguez@uc3m.es
###### Abstract

The translation of a bubble under the action of an acoustic forcing finds applications in fields ranging from drug delivery to sonoluminesce. This phenomenon has been widely studied for cases where the amplitude of the forcing remains constant over time. However, in many practical applications, the duration of the forcing is not long enough for the bubble to attain a constant translational velocity, mainly due to the effect of the history force. Here, we develop a formulation, valid in the limit of very viscous flow and small-amplitude acoustic forcing, that allows us to describe the transient dynamics of bubbles driven by acoustic pulses consisting of a finite number of cycles. We also present an asymptotic solution to this theory for the case of a finite-duration sinusoidal pressure pulse. This solution takes into account both the history integral term and the transient period that the bubble needs to achieve steady radial oscillations, being the former dominant during most of the acceleration process. Moreover, by introducing some additional assumptions, we derive a simplified formula that describes fairly well the time evolution of the bubble velocity. Using this solution we show that the convergence to the steady translational velocity, given by the so-called Bjerknes force, occurs rather slowly, namely as , which explains the slow convergence of the bubble velocity and stresses the importance of taking into account the history force.

Key words:

## 1 Introduction

In 1906, Bjerknes wrote “A pulsating body in a synchronously oscillating current is subject to the action of a resultant force, the direction of which is that of the acceleration in the current at the time when the pulsating body has its maximum volume.” (Bjerknes 1906). The force described here, commonly known as primary Bjerknes force, results from the coupling between the volume oscillations of a body (a bubble in our case) and the pulsation of the external flow. Although the translational motion of a bubble in these circumstances is mainly oscillatory, over many oscillation cycles a net displacement exists, being this the manifestation of the force described by Bjerknes.

Nowadays, there exist different scientific problems and technological applications where the possibility of transporting bubbles via acoustic waves is of great relevance. One of the most promising usages of this force concerns the possibility of using ultrasound to direct microbubbles inside the body towards targets of interest with the purpose of carrying drugs, binding them to specific tissues to aid medical imaging techniques, or to exploit locally the effects of cavitation by making bubbles collapse at precise locations.

In sonochemistry, the appearance of bubbles is usually desirable to promote mixing and to produce localized high temperatures through cavitation. These bubbles have been observed to move in complicated ways in response to sound waves (Rensen et al. 2001), in part due to the not-so-well understood effect of the history force (Toegel et al. 2006).

Another modern technology where the acoustic transport of bubbles could be of interest is the production of hydrogen through artificial photosynthesis (Modestino et al. 2016). In fact, one of the major problems of this technology is the so-called bubble control, that is, how to transport hydrogen bubbles away from the photo-active area once they reach a certain size. One interesting possibility concerns using short acoustic pulses to drive bubbles to desired harvesting locations.

As a last example, and although its implications are still a matter of controversy, we would like to mention the growth of gas bubbles in magna caused by seismic waves (Manga & Brodsky 2006). Indeed, we think that the topic is suggestive and though-provoking enough to deserve a comment here. It the late 1990’s, Sturtevant et al. (1996) and Brodsky et al. (1998) suggested that the triggering of volcanic eruptions by distant earthquakes may be caused by rectified diffusion of gas bubbles in magna caused by traveling seismic waves. Although they took into account rectified diffusion, they did not considered the effect of gas accumulation by coalescence caused by bubble-bubble attraction or bubble-wall interaction. Properly accounting for these effects requires modeling bubble translation in response to short pressure waves.

Motivated by these applications, several authors have studied the problem of the translation of a bubble driven by pulsating pressure waves. A first approach to the problem is phenomenological: that is, to write down a set of equations for the bubble translational dynamics assuming that they are driven by a balance between different effects: drag, inertia of the surrounding liquid (added mass), and Bjerknes force; each of them modeled using an ad-hoc expression. For instance, Rensen et al. (2001) put together a model of this kind to explain q the spiraling behaviour observed in bubbles immersed in an expanding pipe when subject to an acoustic forcing. Dayton et al. (2002) used a similar strategy to describe the motion of a coated bubble in a small pipe, thus the effect of nearby walls must be retained.

A more formal approach was followed by Toilliez & Szeri (2008) who used the method of multiple scales to separate the component of the bubble velocity pulsating at the forcing frequency from the slow drift that yields the net contribution to the mean bubble velocity. Although the authors started their derivation from the unsteady equations of motion for a compressible sphere in the limit of large Reynolds numbers (see Magnaudet & Legendre 1998), they restricted their analysis to the case where the bubble oscillates in a permanent way, once any transient behaviour has damped out.

However, in many practical situations, the duration of the acoustic forcing is not sufficiently large to achieve steady conditions, i.e. a steady mean velocity. In particular, in the case of coated bubbles used in medical applications, short insonation pulses (of the order of 10) are desirable to improve the bubble durability and long-term stability (Chen et al. 2003). In bubble control in artificial photosynthesis, reducing the duration of the pulse while optimizing bubble translation is essential to avoid spending too much of the generated energy in producing the vibration. In these examples, transient effects must be taken into account during the whole duration of the pulse and the subsequent deceleration stage. As we will show in this work, this duration may be as long as several tens or even hundreds of times the viscous time based on the bubble size. Consequently, among other effects, the memory integral appearing in the bubble translation equations (Magnaudet & Legendre 1998) must be accounted for.

In fact, there are several works in the literature that point out the importance of memory effects in the acoustically-driven motion of bubbles. For instance, Toegel et al. (2006) proved that the inclusion of this term was essential to explain the non-periodic, three-dimensional orbits of a bubble around its theoretical equilibrium position observed in single-bubble sonoluminesce experiments when the viscosity of the liquid was increased. In the context of the bubble contrast agents used in medical imaging, Garbin et al. (2009) reported that the trajectory of a microbubble driven by secondary Bjerknes force (the one exerted between two nearby bubbles) could not be properly described unless the history term was included in the dynamic equations.

With these ideas in mind, the aim of this paper is to present a consistent treatment of the acoustically-driven bubble dynamics taking into account transient effects. More specifically, we start from the formulation developed by Magnaudet & Legendre (1998), based on that by Maxey & Riley (1983) but with the viscous terms as given by Yang & Leal (1991), which are valid for a clean gas bubble in the limit of very viscous flow. Then, in order to derive simplified expressions for the evolution of the bubble velocity, we adopt the assumptions that the amplitude of the volume oscillations is small and their period short when measured with the bubble’s viscous time, namely , where is the initial bubble radius and the kinematic viscosity of the liquid.

Although, strictly speaking, the methodology followed here is only applicable in the limit of very viscous flows, we present some experimental evidences suggesting that it is more accurate to retain the memory term and assume zero Reynolds number rather than to neglect this term and model the viscous drag with a constant Stokes’ friction coefficient even if it retains finite-Reynolds effects. This extends the implications of the study to the applications mentioned at the beginning of this section. In fact, although some of them (such as medical imaging) may lay outside the limits of validity of the theory, still the main conclusion holds: the memory force should not be neglected in realistic acoustic pulses.

The paper is thus structured as follows: section 2 shows some experimental evidences that point out the importance of properly accounting for memory effects in the dynamics of bubble translation. Section 3 is devoted to the statement of the problem, and in particular to the presentation of the bubble dynamic equation that will be used in this work. The algorithm used to integrate this equation numerically is then presented in section 4, whereas the asymptotic solution is derived in section 5. Some representative results, obtained through the methodology described previously, are discussed in section 6. Finally, some conclusions are sketched in section 7.

## 2 Experimental evidences

Before starting with the formulation of the problem, it is illustrative to introduce an experimental observation that points out the need to account for the memory effects in the translation of bubbles driven by short acoustic pulses.

We generate isolated hydrogen bubbles with radii, , in the range between 15-25 m by water electrolysis. The cathode is placed at the bottom of a transparent acrylic tank with dimensions 100 40 40 cm. When the bubbles are far enough from the electrode, an ultrasound cylindrical 1-inch transducer with a central frequency of 1 MHz (Sonatest PIM1.0), powered by a pulser-receiver (Ritec RPR-4000), sends a single acoustic pulse of 100 cycles horizontally towards the rising bubble. The usage of a single pulse precludes the occurrence of steady streaming (Riley 2001). This has been checked by adding a low concentration of small tracer particles in some selected experiments. A sketch of the setup is shown in figure 1(a), and more details about the experimental technique can be found in Medina-Palomo (2015). Upon finalization of the experiment, we place the tip of a hydrophone (Ondacorp HGL-0400 with pre-amplifier AH-2010-025) in the field of view of the camera and repeat the acoustic pulse to measure the pressure wave (figure 1).

In response to the acoustic excitation, the bubble departs from its otherwise vertical trajectory moving along the direction of propagation of the acoustic wave, as shown in figure 1(b). Using digital image processing on a high-speed video of the translating bubble acquired with a Memrecam HX-3 operating at 100,000 fps and synchronized with the pulser), we track the trajectory of its center, which allows us to compute the bubble’s horizontal ( component) velocity. A typical plot of this velocity for a pulse is depicted in figures 2(a) and 2(b). We report in table 1 the conditions of five experiments corresponding to bubbles of different sizes translating under the insonation pulse described above. The mean velocity reported in the table is computed using the bubble locations measured along the duration of the pulse (0.1 ms, which, recording at 100,000 fps, yields 10 frames).

For reasons that will be explained in the discussion §6.2, we focus here on the deceleration of the bubble upon the finalization of the pulse. In figure 2 we compare the bubble velocity measured in one of our experiments with the predictions obtained computing the viscous drag with Stokes’ formula:

 Fdrag=−CdμR0V, (2.0)

where is the liquid dynamic viscosity, the equilibrium bubble radius and the mean translational velocity of the bubble. In particular, the equation to be solved is:

 12(4π3ρR30)dVdt=−CdμR0V. (2.0)

Notice that the inertia of the bubble is not given by the (usually negligible) mass of its gas contents, but rather by the added or virtual mass, which is one half of the mass of liquid displaced by a spherical bubble (Leighton 1994). The comparison has been done for two values of the drag coefficient, , namely 4 and 12. The former would correspond to a clean bubble moving steadily at , and it is the smallest value that this coefficient can have. The latter is the value calculated by Magnaudet & Legendre (1998), valid for short times after the deceleration starts. In figure 2(b), where the -axis is plotted in logarithmic scale, it is shown that the experimental values for the deceleration of the bubble cannot be fitted by a straight line. This means that the bubble velocity does not decay exponentially and therefore, the drag coefficient cannot be well described by a constant.

In the light of these results, we can state that the experimental evidence points out that the drag felt by the bubble cannot be modeled by adopting any single value of the Stokes’ drag coefficient. Therefore, we must consider the complete expression of the viscous force, including the history integral. At the end of section 6 we will come back to the analysis of the experiments described here.

## 3 Problem formulation

We consider here the motion of a clean gas bubble of equilibrium radius immersed in a flow that, far away from the bubble, is given by the velocity and pressure fields, and respectively. In particular, we focus on the case where this far field corresponds to an acoustic pressure wave, thus velocity and pressure are related by

 ∂v∞∂t=−1ρ∇p∞. (3.0)

Notice that, although strictly speaking the velocity field given by equation (3) is not incompressible, the fact that it changes along distances of the order of the wavelength, , makes it possible to replace the real compressible flow field by an unsteady but uniform outer flow around the bubble, provided its radius is much smaller than the wavelength, . This assumption, which will be adopted throughout this work, is equivalent to equation (23) of Maxey & Riley (1983), and effectively implies that the outer flow changes over distances much longer than the size of the particle. Additionally, as will be shown below, the smallness of the amplitude of the pressure forcing would also make possible to neglect the nonlinear terms, as they are quadratic in a small parameter.

The starting point for our derivation is the expression given by Magnaudet & Legendre (1998) for the force exerted by the uniform flow of a liquid, with velocity , on a bubble with a time-varying radius that moves with velocity inside the flow. In fact, for a bubble with negligible mass this force must vanish. Consequently, assuming that the Reynolds numbers based on both the bubble’s relative velocity with respect to the fluid and on the wall speed are small, and respectively, we have

 0 = (3.0) +8πνρt∫0exp⎡⎢⎣9νt∫¯tR(s)−2ds⎤⎥⎦erfc⎡⎢ ⎢⎣  ⎷9νt∫¯tR(s)−2ds⎤⎥ ⎥⎦× ×d¯td¯t[R(¯t)(v∞(¯t)−V(¯t))]d¯t,

where the total viscous contribution is given by the last two terms, namely the quasi-steady drag force and the memory integral force. Hereafter, to simplify the notation, the time dependence will be omitted except inside the memory integral term. Notice also that the liquid velocity used in this equation is that particularized at the location of the bubble, so it does not depend on the location anymore. For this reason, we write instead of .

Besides an initial condition for the velocity, to compute the translational response of the bubble, equation (3.0) must be completed with an equation relating the bubble’s volume with the pressure, such as the Keller-Miksis equation (Keller & Miksis 1980):

 (1−˙Rc)R¨R+32˙R2(1−˙R3c)=(1+˙Rc)1ρ(pg−p0)+R˙pgρc−4μρ(˙RR+¨Rc)−2σρR−(1+˙Rc)1ρpac(ω(t−n⋅x/c)) (3.0)

where is the gas pressure inside the bubble, given by

 pg=(p0+2σR0−pv)(RR0)−3κ. (3.0)

Here, the speed of sound in the liquid, the surface tension, the polytropic exponent of the gas evolution inside the bubble, the ambient pressure far away from the bubble in the absence of acoustic waves, and the vapour pressure inside the bubble, which will be neglected henceforth. Finally, the acoustic pressure at the bubble surface is represented by . Notice that we have already introduced here the assumption that the pressure corresponds to that of a traveling wave propagating along the direction given by the unit vector , by letting it depend on . It should be pointed out that, in virtue of the hypothesis , a phase-lag is neglected in . Nonetheless, in the rest of Equation (3.0), terms of order are retained, as they may have a significant contribution in damping bubble volume oscillations (Prosperetti 1977).

### 3.1 Dimensionless formulation

Keeping in mind the motivation of our work, namely the effect of the memory integral term on the translation of bubbles in response to a short acoustic wave, we use the viscous time of the problem, , to define a dimensionless time variable . It will become clear later that the choice of this time scale will simplify the evaluation of the memory integral. The equilibrium radius of the bubble, , is used as length scale, so and . Consistently with this formulation, the pressure field can be expressed as

 pac=pAf({St}τ−Mn⋅~x), (3.0)

with being the amplitude of the pressure wave, the Stokes number and the Mach number. Notice that, under the assumptions adopted here, , since . It is also worth mentioning that the Stokes number that emerges in the problem can be interpreted as the ratio between the viscous time and the period of the acoustic wave.

Introducing equation (3.0) into (3) and making use of the characteristic scales defined above, the corresponding dimensionless equation for the external flow reads

 d~v∞dτ=εM{St}nddτf({St}τ−Mn⋅~x), (3.0)

where the parameter modulates the amplitude of the bubble oscillations. Notice that the spatial derivative found in (3) has been transformed into a temporal one, which will come in handy below in order to integrate the equation.

Similarly, the dimensionless momentum equation yields

 0=εM{St}a3ndf(ξ)dτ+12d(a3~U)dτ+13a~U+23∫τ0G(τ,¯τ)d(a~U)d¯τd¯τ, (3.0)

where we have defined the history kernel as:

 (3.0)

and the relative velocity is hereafter denoted by .

Additionally, to complete the problem, the dimensionless Keller-Miksis equation is required

 ¨a(49M{St}−M{St% }a˙a+a)=12M{St}˙a3+M{St}˙a[(Π0+1{We})a−3κ−Π0]−32˙a2 −49˙aa−{St}2{We}1a+(Π0+1{We}){St}2a−3κ−Π0{St}2−3M{St}κ˙a(Π0+1{We})a−3κ −(1+M{St}˙a){St}2εf(ξ) (3.0)

with being the dimensionless ambient pressure (Euler number) and the Weber number.

As a last step, we rescale the time in order to simplify the treatment of the history integral (Magnaudet & Legendre 1998). This simplification will come in particularly handy for the implementation of the numerical method. Indeed, defining a rescaled time such that, , the kernel (3.0) results

 G(τ,¯τ)=G∗(τ∗−¯τ∗)=exp[τ∗−¯τ∗]erfc[√τ∗−¯τ∗]. (3.0)

Hereafter, this new time variable will be referred to as non-linear time. After some algebra, equations (3.0), (3.0) and (3.0) read

 d~v∞dτ∗=εM{St}ndτ∗dτ∗f({St}τ−Mn⋅~x) (3.0)
 0=εM{St}andf(ξ)dτ∗+d(lna)dτ∗a~U+12d(a~U)dτ∗+13a~U+23τ∗∫0G∗(τ∗−¯τ∗)d(a~U)d¯τ∗d¯τ∗ (3.0)
 ¨a∗(49M{St}−M%St˙a∗a+a)=2(˙a∗)2a(49M{St}−M{St}˙a∗a+a)+12M{St}(˙a∗)3a2 +M{St}˙a∗a2[(Π0+1{We})a−3κ−Π0]−32(˙a∗)2−49a˙a∗−{St}2{We}a3 −Π0{St}2a4+(Π0+1{% We}){St}2a−3κ+4−3M{St}κa2˙a∗(Π0+1{We})a−3κ −(1+M{St}˙a∗a2){St}2a4εf(ξ), (3.0)

where and denote the first and the second derivatives of with respect to the rescaled time .

### 3.2 Some remarks on the underlying hypotheses

We conclude this section with a few comments about the hypotheses adopted thus far and some of their implications. Firstly, as is customarily done in acoustics, the Mach number is assumed to be very small, , which here implies that the size of the particle is much smaller than the length scale of the outer flow, as . Following Maxey & Riley (1983), this allows us to neglect the Faxen correction terms in Equation (3.0) and to approximate in (3).

Secondly, we require the Reynolds numbers based on both the relative velocity between the bubble and the liquid, and on the bubble’s wall velocity, to be small. The smallness of the former justifies the expressions used for the viscous drag (both the steady and history-dependent components), a condition also employed by Maxey & Riley (1983). To express these conditions in terms of the dimensionless parameters defined here, , , and St, we can use order of magnitude analysis on equation (3) to get . Assuming that the bubble’s relative velocity, , is as much of this order, the condition of small translational Reynolds number reads

 R0|v∞−V|ν∼εM{St}≪1. (3.0)

We consider now the condition of small wall-based Reynolds number. The order of magnitude of the amplitude of the radial oscillations is going to denoted by . When the insonation frequency is far from resonance, , whereas when the bubble is excited close to its natural frequency we will have (see Equation (5.0) below). Using this parameter, the condition of small wall-based Reynolds number may be written:

 R0˙Rν∼δ{St}≪1, (3.0)

which translates into in the former case and in the latter. Notice that, close to resonance, this condition results very restrictive at first sight, as it strongly limits the amplitude of pressure wave. Nonetheless, numerical simulations that will be reported elsewhere (Rodríguez-Rodríguez & Martínez-Bazán 2017) suggest that the expressions used here can be applied at least for wall-based Reynolds numbers of order unity. The reason is that the motion induced by the pulsation of the wall is radial and therefore does not generate vorticity. Consequently, although this motion is responsible for part of the virtual mass, its influence on the viscous drag is very limited, at least for order-unity wall-based Reynolds numbers.

Finally, a third restriction imposed by Maxey & Riley (1983) in order to neglect advection effects in the flow around the particle is that , which, using the notation defined here, turns into

 R20|v∞|νλ∼εM2%St≪1. (3.0)

Notice that, effectively, this condition does not imply any further restriction beyond (3.0) and .

Furthermore, although not needed to formulate the problem, the asymptotic approach derived in section 5 will additionally require the bubble radial oscillations to be of small amplitude, , and that their frequency is fast compared to the inverse of the viscous time. In other words, that the Stokes number is large, .

## 4 Numerical method

The translational response of the bubble due to the acoustic pulse is given by equations (3.0), (3.0) and (3.0). They constitute a system of coupled, non-linear, ordinary integro-differential equations that must be solved numerically.

The history integral is computed with the numerical method proposed by Kim et al. (1998), which provides an accurate approximation even when the kernel becomes singular at one of the limits (as is the case for the solid-sphere kernel). More specifically, at the -th time step, the integral is approximated by:

 τ∗∫0G∗(τ∗−¯τ∗)~W¯τ∗d¯τ∗=nΔτ∗∫0G∗(τ∗−¯τ∗)~W¯τ∗d¯τ∗≈ ≈Δτ∗6n−1∑i=1[ G∗(nΔτ∗−(i−1)Δτ∗)~Wτ∗∣∣i−1+ (4.0) +2G∗(nΔτ∗−(i−0.5)Δτ∗)(~Wτ∗∣∣i−1+~Wτ∗∣∣i)+ +G∗(nΔτ∗−iΔτ∗)~Wτ∗|i]+ +0.1Δτ∗2[8√23G∗(0.05Δτ∗)~Wτ∗∣∣n−43G∗(0.1Δτ∗)~Wτ∗∣∣n]+ +0.9Δτ∗6[G∗(Δτ∗)~Wτ∗∣∣n−1+2G∗(0.55Δτ∗)(~Wτ∗∣∣n−1+~Wτ∗∣∣n)+ +G∗(0.1Δτ∗)~Wτ∗∣∣n].

Here, has been denoted by for simplicity. Notice that this method requires the usage of a constant time step, .

Consistently with this, the time marching has been performed with a fourth-order Adams-Bashforth algorithm, which is a multi-step explicit method. The multi-step character of the method facilitates the computation of the memory term, since the algorithm used for its evaluation (Equation 4.0) uses the current, , and previous values. More specifically, the algorithm yields:

 yn+1 = yn+Δτ∗˙y∗(τ∗n,yn) (4.0) yn+2 = yn+1+Δτ∗(32˙y∗(τ∗n+1,yn+1)−12˙y∗(τ∗n,yn)) (4.0) yn+3 = (4.0) yn+4 = (4.0) +3724˙y∗(τ∗n+1,yn+1)−38˙y∗(τ∗n,yn)).

Here, is the vector of unknowns. In the simulations reported in this work, the time step chosen is in the range of . In all of the simulations presented here, we have tested for the temporal convergence of the numerical integration by reducing the time step up to times.

## 5 Asymptotic solutions to the bubble translation problem

As discussed in the previous section, including the history force in the computation of the bubble velocity leads to an integro-differential equation with an integral term that must be evaluated numerically at each time step. Thus, the calculation of the bubble response under a typical acoustic pulse, which might require to integrate for a large number of acoustic cycles, results computationally expensive. With this motivation, we develop in this section an analytical solution for the time evolution of the bubble velocity. Even more importantly, the analytical solution will allow us to obtain physical insight on the physics of the problem and to elucidate which terms are dominant in the bubble dynamics.

We start by assuming that the bubble undergoes small-amplitude volume oscillations, which requires . As mentioned in subsection 3.2, this translates into when the insonation takes place close to the resonance frequency and into otherwise. This poses a difficulty, as the small parameter used in the asymptotic expansions depends on other parameters of the problem (the Stokes number and the dimensionless resonance frequency, , defined below). To avoid this difficulty, and to derive a single formulation that can be used uniformly regardless of the insonation frequency, we expand the bubble radius as

 a=1+εr, (5.0)

although it means that the dimensionless radius, , could be as large as St close to resonance. Notice that this choice does not reduce the validity of the asymptotic analyses performed below, as long as the condition is fulfilled.

Introducing this expression into the dimensionless Keller-Miksis equation (3.0), and neglecting second-order and higher terms, yields:

 (1+4M9{St})d2rdτ2+(49+M{St}Ω20)drdτ+{St}2Ω20r=−{St}2f({% St}τ−Mn⋅~x), (5.0)

with . Furthermore, we can expand the relative velocity times the radius, , assuming that its order of magnitude is dictated by that of , which can be obtained by integrating equation (3.0):

 ~v∞=εM{St}nf({St}τ−Mn⋅~x). (5.0)

Thus, we assume that , which suggests the following ansatz:

 ~W=εM{St}(~W(1)+ε~W(2))+O(ε3). (5.0)

Henceforth, and as has been done when integrating (3.0), we will neglect the influence of the position on the phase of the acoustic wave seen by the bubble, namely the contribution . It can be checked that this phase shift is negligible as long as the bubble translates a distance much smaller than over one acoustic cycle. Furthermore, we will also drop the vector character of the velocities, as they all are aligned with the direction of propagation of the acoustic pulse, .

The smallness of the amplitude of the oscillations also suggests that the difference between the nonlinear and the linear times is going to be small. In order to take advantage of this, we will rewrite Equation (3.0) in terms of and divide the resulting equation by ,

Introducing the expansions for and , and collecting terms of the same order, we get for the leading order, :

 0=dfdτ+L[~W(1)], (5.0)

whereas for the first-order correction, :

 0=rdfdτ+(drdτ−23r)~W(1)+L[~W(2)]. (5.0)

In these equations we have made use of the following linear, but non-local in time, operator:

 L[⋅]=12d[⋅]dτ+13[⋅]+23τ∫0G(τ−¯τ)d[⋅]d¯τd¯τ. (5.0)

It should be pointed out that a term proportional to times the history integral evaluated for has been neglected in (5.0), as its contribution is expected to be small, as will be verified later.

In the following subsections we will derive asymptotic solutions for Equations (5.0), (5.0) and (5.0) for the case of a purely sinusoidal acoustic forcing. In subsection 5.1 we present an approximate solution for an infinitely long acoustic pulse. Besides being directly applicable to describe the acceleration of a bubble in practical situations, this simplified solution will serve us to gain physical insight about which terms of the dynamic equation (5.0) are dominant and thus must be retained to describe the transient motion of the bubble. Then, in subsection 5.2 we lay down the procedure to obtain a solution for the case of a finite-duration pulse of the form

 f(τ)=sin({St}τ)[H(τ)−H(τ−τp)]. (5.0)

Here, is the Heaviside function and the duration of the pulse, which in this study has been chosen to be an integer number of times the period of the acoustic period. Moreover, to further simplify the solution in a limit of relevance for practical applications, the Stokes number is assumed to be large, . That is, the viscous time is much larger than the period of the acoustic wave.

A comment is now in place about the form of the acoustic pulse used in this study (equation (5.0)). Although its amplitude is constant during the duration of the pulse, , for the sake of algebraic simplicity, the theory developed here can be modified to take into account tapering. Indeed, the solution procedure sketched in subsection 5.2 allows to deal with any forcing function that can be decomposed into piecewise linear combinations of trigonometric functions. Thus, it could be straightforwardly applied to pulses tapered at the start with functions of the kind for , with the duration of the tapering, and an analogous expression for the pulse end. We must point out that although the procedure is easy from the conceptual point of view, as it only involves including more segments of the same form (sinusoidal functions) to the definition of the forcing function, the computations would become very tedious.

### 5.1 Simplified solution for the mean velocity in response to a permanent acoustic forcing

As a first assumption, we will consider that the bubble is immediately set to oscillate in a purely periodic fashion when the acoustic pulse starts. This is a reasonable simplification, as the transient component in the solution of Equation (5.0) decays exponentially whereas, as we will show below, the transient terms in the bubble acceleration decay algebraically. Thus, we will approximate the amplitude of the radial oscillations by the permanent solution to the linearized Keller-Miksis equation:

 r=ε{St} ⎷(1+4M9{% St}−Ω20)2{St}2+β2rcos({St}τ+ϕ), (5.0)

where and .

Similarly, to obtain –which is needed to evaluate the second forcing term in Equation (5.0)– we will also neglect the terms arising from the viscous drag in the solution of Equation (5.0), including the history integral. This is justified since they are St times smaller than the time derivatives that appear in the equation. Doing this, we get for the leading-order relative velocity:

 ~W(1)=−2sin({St}τ). (5.0)

Although we do not show it here for the sake of conciseness, it can be checked that this expression agrees very well with the full solution after a few oscillation cycles.

Introducing these expressions into Equation (5.0), and after some algebra, the forcing term may be written

 rdfdτ+(drdτ−23r)~W(1)={St}2(F+Fccos2{St}τ+Fssin2{St}τ) ⎷(1+4M9{St}−Ω20)2{St}% 2+β2r, (5.0)

with

 F = 32cosϕ, (5.0) Fc = −12cosϕ, (5.0) Fs = −12sinϕ. (5.0)

These expressions reveal that the forcing can be decomposed into a constant part plus an oscillation of frequency . As a last simplification, and keeping in mind that we are seeking for the non-oscillatory component of the bubble velocity, we will only consider the constant part of the forcing. Applying the Laplace transform to Equation (5.0), and denoting the amplitude of the constant forcing by

 AF={St}2F ⎷(1+4M9{St% }−Ω20)2{St}2+β2r, (5.0)

we may write:

 0=AFs+(12s+13+23s√s(1+√s))~\@fontswitchW(2), (5.0)

where is the Laplace transform of . Equation (5.0) can be inverted exactly, as pointed out by Stepanyants & Yeoh (2009) for an analogous case, yielding for the relative velocity:

 ~W(2)=AF{Fs1+e−ατ(Ds1cos(γτ)+Es1−Ds1αγsin(γτ))+Cs1es1τerfc√s1τ (5.0)

Here, , , and are numerical constants that have been defined in Appendix B, whereas functions and are specified in Appendix A.

Finally, we recall the definition of the radius-weighted relative velocity, , to obtain the bubble mean velocity. Indeed, keeping in mind that has zero mean, and assuming , we get for the non-oscillatory component of the velocity:

 ~Vno≈−ε2M{St}~W(2). (5.0)

Notice that approximating the radius by unity does not modify the non-oscillatory part of the velocity at leading order, as is a purely oscillatory function with the approximations made in this subsection. Expressions (5.0) and (5.0) are compared with the numerical solution of the full problem in figures 3, 4 and 5. Also, the time evolution of the bubble radius and its time derivative is plotted in figure 6. The discussion of the behavior of this solution is postponed to §6, to continue now with the derivation of a more complete asymptotic solution of the problem.

### 5.2 Response to a pulse of finite duration

Here we consider the response to an acoustic pulse consisting in a sinusoidal forcing of frequency starting at and ending at a finite time , which corresponds to an integer number of oscillation periods, as defined in Equation (5.0). For the sake of conciseness, and due to the length of the full solution, we have relegated its expression to Appendixes D, E and F. In this subsection we just point out the method followed in its derivation.

We start with the solution to the leading order component, namely , which is mainly oscillatory. Applying Laplace transform to Equation (5.0), yields

 0=\@fontswitchF(s)+12(s~\@fontswitchW(1)(s)−~W(1)|0)+13~\@fontswitchW(1)(s)+23⎛⎜⎝s~\@fontswitchW(1)(s)√s(1+√s)⎞⎟⎠, (5.0)

with the Laplace transform of , the Laplace transform of the leading order of the radius-weighted relative velocity, , and its initial value. Isolating the Laplace transform of the radius-weighted relative velocity:

 ~\@fontswitchW(1)(s)=~W(1)|0(s2+s−23+43√s)(s−s1)[(s+α)2+γ2]+\@fontswitchF(s)−2(s2+s−23+43√s)(s−s1)[(s+α)2+γ2]. (5.0)

The solution has been split into two components. The former, , is due to the initial relative velocity, non-zero in case of having a bubble immersed in a flow that is not at rest before the acoustic forcing begins, whereas the latter, , is due to the acoustic forcing and thus depends on the function describing the shape of the acoustic pulse.

To keep advancing, it is necessary to specify the forcing function . Recalling Equation (5.0), its derivative reads

 df(τ)dτ={St}cos({St}τ)[H(τ)−H(τ−τp)]+sin({St}τ)[δ(τ)−δ(τ−τp)]. (5.0)

Since the duration of the pulse, , is chosen here to be an integer number of times the period of the acoustic wave, the term containing Dirac’s deltas is zero. Therefore, will be given by

 df(τ)dτ={St}cos({St}τ)[H(τ)−H(τ−τp)] (5.0)

and thus

 (5.0)

Substituting this expression into (5.0), and inverting the Laplace transforms we obtain , whose expression is provided in Appendix D.

We focus now on the non-oscillatory part of the relative velocity, which mainly arises from the first-order correction, , given by Equation (5.0). Indeed, this component is the responsible for the appearance of a net force, namely the Bjerknes force. Notice that this equation has two forcing terms: one due to the acoustic forcing and another due to the contribution of the leading-order component, . As the equation is linear in , we can search for the particular solutions associated to both contributions separately:

 0=rdfdτ+L[~W(2)1] (5.0)
 0=(drdτ−23r)~W(1)+L[~W(2)2] (5.0)

with .

In the case of equation (5.0), the forcing term is a sum of products of decaying exponentials and trigonometric functions, modulated in time by the Heaviside step function. Therefore, can be obtained by applying the Laplace transform method as explained in subsection 5.1.

However, in the case of equation (5.0), the forcing term is not so simple. Apart from the trigonometric and exponential terms, the Faddeeva function as well as functions and also appear in it. Consequently, the Laplace transform method cannot be directly applied, as the inverse Laplace transform arising from this component cannot be expressed in terms of known functions. To circumvent this difficulty, we replace the Faddeeva, Isin and Icos functions by suitable approximations derived in Appendix A. More exactly, only the terms representing a permanent oscillation have been retained. The complete analytical solutions of components and are given in appendixes E and F respectively.

Once the radius-weighted relative velocity is known, the actual bubble velocity is easily obtained, since the analytical solution of the external flow is also known. Indeed, for the forcing considered,

 ~v∞(τ)=εM{St}sin({St}τ)[H(τ)−H(τ−τp)], (5.0)

so the bubble velocity finally results:

 ~V(τ)=εM{St}sin({St}τ)[H(τ)−H(τ−τp)]−εM{St}(~W(1)(τ)+ε~W(2)(τ))a(τ). (5.0)

## 6 Results and Discussion

### 6.1 Numerical results

We compare here the bubble velocities predicted by the numerical solution of the full problem formulated in §3 with those obtained from a simplified version of the system without the history integral. Moreover, we also employ the asymptotic approach presented in the previous section to evaluate its accuracy versus the full solution. For the sake of conciseness, three representative cases have been selected, whose conditions are reported in table 2. They are all close to resonance, since it is when the bubble moves at the highest speed, which is the situation desired in practical applications. We would like to stress here that even though the number of cycles is relatively large in these cases (), still memory effects are important to compute translation velocities. Were we using shorter pulses, like those used in medical imaging (), these effects would be even more important.

The results are plotted in figures 3, 4 and 5. Panels labelled as () show the oscillatory velocities corresponding to the numerical solution of the problem (blue solid lines) and to the asymptotic solution (cyan dashed lines). More interestingly, panels labelled as () contain the period-averaged velocities, defined as

 ⟨~V⟩={St}2π∫τ+2π/{% St}τ~V(s)ds. (6.0)

In these panels, the period-averaged velocities obtained from the numerical solution of the problem without taking into account the history integral are also represented with black dot-dashed lines. As the experimental observations described in section 2 suggest, our calculations confirm that the proper description of the bubble velocity requires taking into account the memory integral, even though the pulses last for several tens of viscous times, i.e. a time much longer than the viscous one. Indeed, figures 3(), 4(), and 5() show that the period-averaged velocity calculated neglecting history achieves the terminal value after about ten viscous times, a period much shorter than the duration of the pulse. On the contrary, the velocity predicted by the full solution approaches the asymptote much more slowly, not even reaching it while the acoustic excitation is on. In figures 3(), 4(), and 5() the solutions have also been extended beyond the duration of the acoustic pulse, to illustrate the slow convergence towards the terminal velocity. A similar trend is observed once the pulse stops: the velocity does not decay exponentially, as would be expected if no memory integral was considered, but rather potentially.

For the sake of completeness, the first few cycles of the radial oscillations are also shown in figure 6 for the three cases considered, along with their time derivative. The behavior of the full numerical solution for the bubble velocity is fairly well captured by the asymptotic solution derived in subsection 5.2, specially that of the period-averaged velocity, which is the responsible for the net bubble translation. This is quite remarkable taking into account that the cases reported in this work are all relatively close to the resonance frequency of the bubble. We have chosen cases close to resonance since they imply that the amplitude of the radial oscillations is going to be largest, which translates into strongest Bjerknes force, and thus to very large bubble displacements. In other words, this means that it is near resonance when the formulation developed here is going to be of greatest interest for the applications. Although near resonance the amplitude of the radial oscillations may not be small in general, we impose an additional requirement to the parameters of the problem, namely , which translates into when the oscillations occur close to resonance (see Equation (5.0)). It is worth paying attention to case II, where this condition is nearest to not being satisfied (). It is interesting to observe how, even if the amplitude of the oscillations is not properly captured, still the average velocity is barely affected. This is consistent with the fact that, in the weakly non-linear regime, the response of the bubble differs from the linear solution in the appearance of harmonics and sub-harmonics that do not contribute to the average of the forcing terms in Equation (5.0).

The simplified solution presented in subsection 5.1 (short-dashed red line in figures 3(), 4(), and 5()) also describes fairly well the averaged solution of the full problem while the acoustic forcing is on. In particular, once the transient behavior in the bubble radial oscillations has damped out. Thus, although the exact evaluation of the bubble velocity taking into account history effects may be computationally expensive or analytically involved, the procedure followed in subsection 5.1 allows for an accurate description of the bubble dynamics while keeping mathematical or computational complexity within reasonable limits. Interestingly, it seems that to capture the transient bubble dynamics it suffices to retain the memory integral only at first order, whereas the transient volume oscillations and the history integral at leading order seem to influence the bubble velocity only at very short times.

It is reasonable to expect that this simplified solution may not describe well the dynamics of a bubble insonated with a pulse consisting of a very small number of cycles, in particular when the duration of the pulse is smaller or of the order of a few viscous times, where the volume oscillations have not yet reached the permanent regime. Notice that this limitation only affects the simplified solution described in §5.1, not the more rigorous asymptotic approach presented in §5.2, which properly accounts for the effect of the transient volume oscillations. It should be pointed out that numerical simulations performed for the same conditions reported in table 2 but with a duration of ten pulses reveal that the history force is still very important and that the solution with no history term completely overpredicts the bubble velocity after the very first cycles.

Moreover, the good agreement between the numerical results and the simplified solution suggests that it can be used to infer how the bubble velocity approaches the steady translational velocity. To that end, we introduce into (5.0) the asymptotic approximations of Appendix A, which results into

 ~Vno≈−ε2M{St}AF(Fs1+Cs1/√s1−8Ds2/3√πτ), (6.0)

neglecting terms that decay exponentially. Consequently, we can conclude that the bubble net translational velocity approaches the steady one as rather than exponentially, as would be the case were the memory effects neglected. This explains the slow convergence observed in our numerical simulations and supports the conclusion that history is indeed very important in order to accurately compute the bubble dynamics in many practical situations.

Finally, we would like to point out that most of the conclusions mentioned here concerning the rate at which the velocity converges to an asymptote are of application for bubbles in a fluid contaminated with surfactants. Indeed, in such case, the kernel that must be used in the memory integral is that of a solid sphere, namely . Using the approximations of appendix A, it can be shown that the kernel for a clean bubble converges to this one at long times.

### 6.2 Comparison with experiments

We come back now to the experimental results that motivated this study. We compare the experimental deceleration of the bubble after the pulse (see section 2) with the predictions obtained using the two extreme values of the Stokes’ coefficients ( and ) and with the analytical solution considering the effect of the history force.

We focus here on the deceleration period, as a direct comparison of the bubble dynamics during the pulse is not feasible at the large pressure amplitudes used in the experiments. In fact, the radial dynamics at these pressures are not well predicted by any equation that assumes radial symmetry, as shape oscillations may appear. However, these high pressures are needed in our experiments for the displacements to be large enough to be accurately measured. Thus, we compare instead the asymptotic trend that the bubble velocity follows at long times upon the pulse is off with the experimental results.

It can be easily seen that, in absence of forcing, the time evolution of the bubble velocity is given by Equation (D 0), but replacing by , namely the velocity at the moment the pulse ends:

 ~V(τ)=~V|0{e−ατ(A0cos(γτ)+B0−A0αγsin(γτ))+C0es1τerfc√s1τ++43[(E0−2αD0)Icos(τ,α,γ)+−E0α+D0(α2−γ2)γ