1. Introduction
With their large number of applications in natural phenomena, technology and medicine, gas bubbles in liquids have been the object of a voluminous literature. In spite of their apparent simplicity, these minute entities present unexpectedly complex features that are not easily incorporated into models. In particular, the dynamic response of gas bubbles to variations of the ambient pressure depends on the gas pressure in the bubble, which, in its turn, is critically dependent on the gas temperature. In addition to its effect on the gas pressure, the liquid–gas energy exchange is also the dominant mechanism for the conversion of compression work into heat, thereby determining the dissipation of mechanical energy.
These processes are well understood in the linear approximation (see e.g. Devin Reference Devin1959; Plesset & Prosperetti Reference Plesset and Prosperetti1977). In the more interesting nonlinear case, however, while a reasonable theory has been available for some time (see e.g. Prosperetti Reference Prosperetti1991), it has not been possible so far to simplify it to the point of making it widely usable without a relatively large investment of effort. For this reason, nearly universal use is made in the literature of the polytropic model, in which the gas pressure 
$p$ is approximated in terms of the bubble volume 
$V$ as
\begin{equation} pV^\kappa = p_0V_0^\kappa ,\end{equation}
with 
$\kappa$ a polytropic index and the subscript 0 denoting equilibrium values; in particular, 
$V_0=\frac {4}{3}{\rm \pi} R_0^3$ with 
$R_0$ the equilibrium radius. While this model may represent the real behaviour of the gas pressure in restricted conditions and parameter ranges, it can be used with confidence only within the narrow boundaries of small-amplitude motion with a well-defined frequency and, furthermore, it completely misses the thermal dissipative effects. Another important deficiency is the inability to reliably predict the bubble temperature, thus limiting its usefulness in applications such as sonochemistry. Furthermore, the stability of the bubble and, ultimately, its integrity depend on the radial acceleration, which rests on the correct prediction of the internal pressure. A more sophisticated model was developed by Toegel, Hilgenfeldt & Lohse (Reference Toegel, Hilgenfeldt and Lohse2002) in terms of the average temperature in the bubble. While this model gives a better agreement with the detailed theory (Stricker, Prosperetti & Lohse Reference Stricker, Prosperetti and Lohse2011), it provides no information on the temperature distribution in the bubble, which can be strongly non-uniform.
It is the purpose of this paper to present a simplified version of the theory mentioned before, which, in spite of some limitations, accurately reproduces the correct gas pressure and temperature in a useful range of conditions and parameters. Crucially, this approximate theory is cast in the form of ordinary differential equations, which therefore do not add much of an overhead to the modelling. The simplified theory is compared with a full Navier–Stokes solution of the thermo-fluid dynamics of the gas bubble to demonstrate its accuracy and to identify its limitations.
Since the focus here is on the gas contained in the bubble (except in one case, see figure 5), for the radial motion we use a simple incompressible Rayleigh–Plesset equation, namely
\begin{equation} R\ddot{R}+\frac{3}{2}\dot{R}^2 = \frac{1}{\rho_l}\left[p(t)-P_l(t)- \frac{2\sigma}{R}-4\mu_l\frac{\dot{R}}{R}\right] ,\end{equation}
with 
$\rho _l$ and 
$\mu _l$ the liquid density and dynamic viscosity and 
$\sigma$ the surface tension coefficient. The instantaneous bubble radius is 
$R=R(t)$, the gas pressure in the bubble is 
$p(t)$ and the time-dependent ambient pressure in the liquid is 
$P_l(t)$; the overdot denotes time differentiation.
2. Gas pressure and temperature in the homobaric approximation
A simple estimate based on the integration of the gas momentum equation shows that the pressure difference between the bubble centre and its surface is of the order of
\begin{equation} \frac{p(R(t),t)-p(0,t)}{\bar{p}} = O\left(\frac{R}{\lambda} \frac{\dot{R}}{c},\frac{\dot{R}^2}{c^2}\right) , \end{equation}
with 
$\bar {p}$ an average pressure, and 
$\lambda$ and 
$c$ representative values of the sound wavelength and speed in the bubble. The smallness of the quantities on the right-hand side of this estimate suggests that the spatial variation of 
$p$ in the bubble can be neglected. As shown in Prosperetti (Reference Prosperetti1991), with the assumption of a perfect nature of the gas in the bubble, this observation leads to a differential equation for the gas pressure:
\begin{equation} \dot{p}=\frac{3}{R} [(\gamma -1) k[\partial_rT]_{r=R}-\gamma p \dot{R} ], \end{equation}
in which 
$T$ is the absolute temperature, 
$\gamma$ the ratio of the gas specific heats and 
$k$ the gas thermal conductivity. A second consequence is a differential equation for the gas temperature, which, in terms of the scaled radial coordinate 
$y=r/R(t)$, is
\begin{equation} \frac{\gamma}{\gamma-1}p\partial_tT + \frac{1}{R^2} (k\partial_yT-y[k\partial_yT]_{y=1} ) \partial_yT = T\dot{p} + \frac{T}{R^2y^2}\partial_y (k y^2\partial_y T). \end{equation}The second group of terms on the left-hand side represents the convective contribution to gas energy transport. This partial differential equation is not particularly difficult to solve (see e.g. Kamath & Prosperetti Reference Kamath and Prosperetti1989), but including this step in the modelling of complex bubble phenomena is cumbersome and has not been frequently done.
In basing our simplified models on this formulation, we will keep 
$T_s$, the liquid temperature at the bubble surface, constant. This step is not necessary, as one could couple the energy equation in the liquid to our approximate model, but is in keeping with the aim to develop a simplified theory. As for the accuracy of the approximation, a simple estimate of the magnitude of the variation of 
$T_s$ is given by (see e.g. Kamath & Prosperetti Reference Kamath and Prosperetti1989) 
$(T_s-T_\infty )/(T_c-T_s) \simeq \sqrt {k_gc_{pg}\rho _g/(k_lc_{pl}\rho _l)}$, with 
$T_\infty$ the liquid temperature far from the bubble, 
$T_c$ the gas temperature at the centre of the bubble, 
$c_p$ the specific heat, and indices 
$g$ and 
$l$ referring to gas and liquid, respectively. With typical values of the parameters involved, we find that 
$(T_s-T_\infty )/(T_c-T_s) \sim 10^{-3}$, which justifies the approximation (see also Stricker et al. Reference Stricker, Prosperetti and Lohse2011).
3. Simplified models
An important constraint to be satisfied in developing a simplified theory of thermal processes in the interior of a bubble is conservation of the gas mass: unless mass is strictly conserved, the bubble exhibits a spurious growth or shrinkage (Prosperetti Reference Prosperetti1991).
In the perfect-gas approximation, 
$m_0$, the total mass of gas initially present in the bubble, can be written as 
$m_0=\frac {4}{3}{\rm \pi} R_0^3p_0/({\mathcal {R}}_GT_0/M)$, with 
${\mathcal {R}}_G$ the universal gas constant and 
$M$ the molecular mass of the gas. The invariance of this mass in the course of the bubble motion requires the constancy of
\begin{equation} m_0 = 4{\rm \pi} \int_0^{R(t)} r^2 \rho(r,t) \,\textrm{d} r= \frac{4{\rm \pi} R^3 p}{{\mathcal{R}}_G/M} \int_0^1 \frac{y^2}{T}\,\textrm{d} y ,\end{equation}
in which 
$\rho$ is the gas density. It is not difficult to show that the time derivative of 
$m_0$ as expressed by this relation vanishes exactly if 
$T$ satisfies the energy equation (2.3), which suggests that this equation is a good starting point for the development of an approximate model. We will use it to calculate the gas pressure from a suitable approximation to the temperature distribution.
After multiplication of the energy equation (2.3) by 
$y^2$ and integration over the bubble volume, a few simple steps lead to the result
\begin{align} \frac{\gamma}{\gamma -1} p \frac{\textrm{d}}{\textrm{d} t}\langle T \rangle = -\frac{3\gamma p \dot{R}}{R} \langle T \rangle+\frac{3k}{R^2}\left[ [2T_s +(\gamma-2) \langle T \rangle ][\partial_y T]_{y=1} - 2\int_0^1 y^2 (\partial_yT)^2\, \textrm{d} y \right] ,\end{align}
in which we have used (2.2) to express 
$\dot {p}$ in the first term on the right-hand side and 
$\langle T \rangle$, the volume-averaged gas temperature, is defined by
\begin{equation} \langle T \rangle = 3\int_0^1 y^2 T\, \textrm{d} y .\end{equation}In the derivation we have neglected the temperature dependence of the gas thermal conductivity, for which the value corresponding to the undisturbed liquid temperature should be used, given its importance for the determination of the heat flux at the wall.
We developed two versions of the simplified model. In the first one the gas temperature is approximated as a biquadratic function
\begin{equation} T(y,t) = T_c(t) + A(t) y^2 + B(t)y^4 ,\end{equation}
with 
$T_c$ the centre temperature. Owing to the constancy of the bubble surface temperature, it is necessary that 
$T_c+A+B=T_s$ at each instant, so that the temperature in this model has, in fact, two degrees of freedom. The simpler model is quadratic, and can be obtained from the biquadratic one simply by taking 
$B=0$ and dropping equation (3.6) below. This particular ansatz for the temperature distribution is inspired by the nature of the solutions of the energy equation in the limit of large thermal diffusivity derived in Prosperetti (Reference Prosperetti1991). In that work it is shown that the first correction to an isothermal gas is a quadratic term in 
$y$, while the second correction contains terms proportional to 
$y^2$ and 
$y^4$. It is also interesting to note that a simple quadratic polynomial happens to be an exact solution of the energy equation (2.3). This solution is, however, of little interest because, for consistency of the relations to which it leads, it requires that 
$A = 0$ so that the gas remains isothermal. We have also experimented with a bicubic polynomial of the form 
$T(y,t)=T_c+Cy^3+Dy^6$, which has the advantage of simplifying the integration in (3.1) but otherwise performs similarly to (3.4).
In the end we opted for the biquadratic form, as the numerical method used in several earlier studies for a precise integration of the energy equation is based on an expansion of the unknown temperature in a series of even Chebyshev polynomials, which only contain even powers of 
$y$ (see e.g. Kamath & Prosperetti Reference Kamath and Prosperetti1989; Hao, Zhang & Prosperetti Reference Hao, Zhang and Prosperetti2017). Additional powers of 
$y$ could be added to (3.4) with several options for the equations necessary to determine the new coefficients. One possibility would be a collocation method, in which the energy equation (2.3) is enforced at a suitable number of collocation points 
$0\leq y_j<1$. Indeed, one version of the Chebyshev expansion method could be considered as an example of this approach. Another possibility would be to use higher-order moments of the temperature distribution beyond (3.3) such as 
$\int _0^1 y^{n+2}T(y,t)\,\textrm {d} y$ with 
$n=1,2,\ldots$. These extensions would, however, render necessary the evaluation of the integral in (3.1) by numerical, rather than analytical, means. Thus, as long as a simple approximate model is acceptable, it seems that the use of the biquadratic approximation recommends itself in terms of its simplicity and, as will be shown, its reasonable accuracy.
In order to determine one of the two degrees of freedom of the temperature distribution (3.4), we use the integral form (3.2) of the energy equation, finding
\begin{align} \frac{\gamma p}{\gamma -1} \frac{\textrm{d}}{\textrm{d} t}\langle T \rangle = -\frac{3\gamma p \dot{R}}{R} \langle T \rangle\!+\!\frac{6k}{R^2}\left[ [2T_s \!+\!(\gamma-2) \langle T \rangle ](A\!+\!2B) -\dfrac{4}{5}A^2-\dfrac{16}{9}B^2 - \dfrac{16}{7} AB \right] ,\end{align}
with 
$\langle T \rangle = T_c+\frac {3}{5}A+\frac {3}{7}B$. To determine the other degree of freedom, we use the differential form of the energy equation evaluated at the bubble centre, which, after eliminating 
$\dot {p}$ by means of (2.2), is
\begin{equation} \frac{\gamma}{\gamma-1} p \dot{T}_c=\frac{3T_c}{R^2} [2(\gamma -1) k (A+2B)-\gamma p R \dot{R}] +\frac{6kT_c A}{R^2}.\end{equation} The pressure is calculated from the expression (3.1) for the gas mass. The calculation of the integral of 
$1/T$ in (3.1) gives the following results, in all of which the expression 
$[\tan ^{-1}\sqrt {X}]/\sqrt {X}$ should be interpreted as 
$[\tanh ^{-1}\sqrt {-X}]/\sqrt {-X}$ for 
$X<0$:
(i) If
$h=A^2-4T_cB>0$, with $2f = A - \sqrt {h}$ and $2g = A+\sqrt {h}$, we find
(3.7)\begin{equation} \int_0^1 \frac{y^2}{T_c+Ay^2+By^4}\, \textrm{d} y = \frac{1}{\sqrt{h}} \left[\sqrt{\frac{g}{B}} \tan^{-1}\sqrt{\frac{B}{g}} - \sqrt{\frac{f}{B}} \tan^{-1}\sqrt{\frac{B}{f}} \right] .\end{equation}(ii) If
$h=A^2-4T_cB<0$, the result is
(3.8)in which the real and imaginary parts of the complex quantity\begin{align} \int_0^1 \frac{y^2}{T_c+Ay^2+By^4}\, \textrm{d} y = \frac{1}{\sqrt{-h}}\left[ \varGamma_i \log \left(\frac{\sqrt{T_c}+\sqrt{B}+\sqrt{2\sqrt{BT_c}-A}}{\sqrt{T_s}} \right) + \theta \varGamma_r \right] , \end{align}$\varGamma =\varGamma _r-\textrm {i}\varGamma _i$ are given by
(3.9a,b)and the angle\begin{equation} \varGamma_r = \sqrt{\frac{1}{2}\left(\sqrt{\frac{T_c}{B}}- \frac{A}{2B}\right)} , \quad \varGamma_i = \sqrt{\frac{1}{2}\left(\sqrt{\frac{T_c}{B}}+ \frac{A}{2B}\right)} , \end{equation}$\theta$ is found from $\sqrt{T_s}\,\cos\theta = \sqrt {T_c}-\sqrt {B}$, $\sqrt{T_s}\,\sin\theta = \sqrt {A+2\sqrt {T_cB}}$.(iii) If
$B=0$ and $A>0$,
(3.10)\begin{equation} \int_0^1 \frac{y^2}{T_c+Ay^2+By^4}\, \textrm{d} y = \frac{1}{A}\left[1- \sqrt{\frac{T_c}{A}}\tan^{-1}\sqrt{\frac{A}{T_c}} \right] .\end{equation}(iv) If
$h=0$ and $A>0$,
(3.11)while, for\begin{equation} \int_0^1 \frac{y^2}{T_c+Ay^2+By^4}\, \textrm{d} y = \frac{1}{2B} \left[ \left(\frac{B}{T_c}\right)^{1/4}\tan^{-1} \left(\frac{B}{T_c}\right)^{1/4} -\frac{1}{1+\sqrt{T_c/B}} \right] ,\end{equation}$A<0$,
(3.12)\begin{equation} \int_0^1 \frac{y^2}{T_c+Ay^2+By^4}\, \textrm{d} y = -\frac{1}{2B} \left[ \left(\frac{B}{T_c}\right)^{1/4}\tanh^{-1} \left(\frac{B}{T_c}\right)^{1/4} +\frac{1}{1-\sqrt{T_c/B}} \right] .\end{equation}
In order to ascertain the validity of the models, we compare their results with those of a finite-volume-based numerical solution of the spherically symmetric balance equations for the gas mass, momentum and total energy coupled with the Rayleigh–Plesset equation (1.2) for the motion of the boundary. The bubble volume is discretized into 
$N$ spherical shells indexed by 
$i=1, 2, \ldots , N$, each one with volume 
$V_i=\frac {4}{3}{\rm \pi} R^3(y^3_{i+1/2}-y^3_{i-1/2})$ bounded by surfaces of area 
$S_{i+1/2}=4{\rm \pi} R^2 y_{i+1/2}^2$, with 
$y_{i+1/2}R$ denoting the surface separating 
$V_i$ from 
$V_{i+1}$; we define 
$y_{1/2} = 0$ and 
$y_{N+1/2} = 1$. The equations are discretized by a finite-volume method. The discretized form of the equation of continuity, for example, is
\begin{equation} \dot{M}_i = F_{i-1/2}-F_{i+1/2} , \end{equation}
with 
$M_i= \rho _i V_i$ the gas mass in the spherical shell. The fluxes are given by 
$F_{i+1/2}= \rho _{i+1/2}^{uw} [u_{i+1/2}- y_{i+1/2} \dot {R}] S_{i+1/2}$, with the density 
$\rho _{i+1/2}^{uw}$ approximated with a second-order-accurate upwind interpolation, which is found to be necessary for numerical stability during the violent phases of the bubble motion. Subtraction of the interface velocity 
$y_{i+1/2} \dot {R}$ is necessary since the use of the variable 
$y$ causes the interfaces to move with the changing bubble radius. Use was made of 50 cells with refinement near the interface; standard convergence studies proved that this resolution was adequate.
Summary of the approximate models: For the convenience of the reader, we summarize here the equations that need to be solved for the present models, both of which rely on a suitable form of the Rayleigh–Plesset equation (compressible or incompressible) to couple the bubble radius with the internal gas pressure.
(a) For the biquadratic model, the two degrees of freedom of the gas temperature are determined by integrating the two ordinary differential equations (3.5) and (3.6). Any two of
$\langle T \rangle$, $T_c$, $A$ or $B$ can be taken as the unknowns, noting the relations $T_c+A+B=T_s$, with $T_s$ the constant liquid temperature, and $\langle T \rangle = T_c+\frac {3}{5}A+\frac {3}{7}B$. The pressure is found from (3.1) in which the integral is to be evaluated according to one of (3.7)–(3.12), whichever is appropriate.(b) For the quadratic model, the gas temperature has only one degree of freedom that can be determined by integrating (3.5) with
$B = 0$ and $A = T_s-T_c$. The pressure is found from (3.1) in which the value of the integral is given by (3.10).
4. Results
In illustrating the performance of the approximate models, we will study bubbles driven into oscillation by a harmonically varying ambient pressure
\begin{equation} P_l(t) = p_\infty - p_a \sin \omega t , \end{equation}
with 
$p_\infty$ the static pressure, taken as 100 kPa, and 
$p_a$ the acoustic pressure amplitude oscillating at the angular frequency 
$\omega$. Since the focus of this work is on the bubble interior, we will not vary the physical properties of the liquid, just using those of water at 
$T_s = 20\,^\circ$C with 
$\rho _l = 998$ kg m
$^{-3}$, 
$\sigma = 0.0725$ J m
$^{-2}$ and 
$\mu _l = 10^{-3}$ Pa s. For the gas we use the properties of air at the same temperature with 
$k = 0.0259$ W m
$^{-1}$ K
$^{-1}$.
A good way to gain a quick overall impression of the models’ performance is a consideration of graphs of the normalized maximum oscillation amplitude,
\begin{equation} R_M^*= \frac{R_{max}-R_0}{R_0} , \end{equation}
with 
$R_{max}$ the maximum radius attained in the course of a steady oscillation, versus the ratio 
$\nu /\nu _0$ with 
$\nu =\omega /2{\rm \pi}$ and 
$\nu _0$ the natural frequency for linear oscillations of the bubble. The values of 
$\nu _0$ for the bubble radii that we consider are given in table 1 together with the linear-theory polytropic index.
Equilibrium bubble radius 
$R_0$, linear resonance frequency 
$\nu _0$ and linear polytropic index 
$\kappa$ for the bubble radii considered in this work.
The general appearance of graphs of 
$R_M^*$ versus frequency, introduced by Lauterborn (Reference Lauterborn1976), is well known. Resonance peaks appear in the neighbourhood of 
$\nu /\nu _0$ equal to the ratio of two small integers corresponding to harmonic, subharmonic and ultra-harmonic resonances. Since the bubble is a ‘softening’ oscillator, the frequency resulting in the largest oscillation amplitude decreases with increasing amplitude and, as a consequence, the resonance peaks bend to the left unlike the case of a linear oscillator. This feature results in frequency ranges in which the steady oscillation amplitude has two stable values, located on the solid portion of the lines sketched in figure 1, and one unstable value located on the dashed portion of the line. Which steady state will ultimately be attained in any specific case depends on the fact that the space of initial conditions is divided into two domains of attraction, one for each stable steady state.
Illustrating the general features of a resonance peak for a nonlinear softening oscillator such as a gas bubble; the dashed lines show unstable points.
In the present study, all the simulations were started from equilibrium. In the direction of increasing frequency, we found that the oscillation amplitude gradually increased following the lower rising branch of the resonant peak of figure 1. This behaviour indicates that the state of equilibrium belongs to the domain of attraction of the lower stable steady state. When the point of vertical tangency of the curve 
$R_M^*$ versus 
$\nu /\nu _0$ was reached, the steady oscillation amplitude abruptly jumped to the other larger branch of the resonant peak, as shown by the upward (blue) arrow in figure 1. As the frequency increased further, the amplitude decreased following this descending branch until the next resonance was approached and a similar behaviour repeated itself. Realizing the jump indicated by the downward (red) arrow of figure 1 would require initial conditions placed in the domain of attraction of the larger stable steady state.
Figure 2 shows 
$R_M^*$ versus 
$\nu /\nu _0$ for different bubble radii. The left and right panels show results obtained with the biquadratic and quadratic temperature distributions, respectively. Panels (a) and (e) are for bubbles with 
$R_0 = 1\ \mathrm {\mu }$m with 
$p_a/p_\infty = 1.6$; panels (b) and (f) for 
$R_0 = 10\ \mathrm {\mu }$m with 
$p_a/p_\infty = 0.8$; panels (c) and (g) for 
$R_0 = 100\ \mathrm {\mu }$m with 
$p_a/p_\infty = 0.5$; and panels (d) and (h) for 
$R_0 = 500\ \mathrm {\mu }$m with 
$p_a/p_\infty = 0.5$. The red dots connected by red lines are the Navier–Stokes results, while the blue symbols are from the approximate models.
Normalized maximum radius 
$R_M^*$ (4.2) for steady oscillations versus sound frequency 
$\nu$ normalized by the bubble natural frequency 
$\nu _0$. Red lines and points are Navier–Stokes simulations; blue symbols are approximate models, biquadratic in the left column and quadratic in the right column. Equilibrium radius 
$R_0 = 1\ \mathrm {\mu }$m and normalized acoustic pressure amplitude 
$p_a/p_\infty = 1.6$ in (a) and (e), 
$R_0 = 10\ \mathrm {\mu }$m and 
$p_a/p_\infty = 0.8$ in (b) and (f), 
$R_0 = 100\ \mathrm {\mu }$m and 
$p_a/p_\infty = 0.5$ in (c) and (g), and 
$R_0 = 500\ \mathrm {\mu }$m and 
$p_a/p_\infty = 0.5$ in (d) and (h).
Starting from the results for 
$R_0 = 1\ \mathrm {\mu }$m in the first row of figure 2, we see that both models capture the bubble behaviour very well. As indicated by the value of the linear polytropic index in table 1, such small bubbles are very nearly isothermal and the temperature field is well approximated by both models. The situation is similar for 
$R_0 = 10\ \mathrm {\mu }$m in the second row. In panel (f), the isolated point near the peak of the first harmonic corresponds to an ultra-harmonic oscillation of order 6/3. Recovering the harmonic Navier–Stokes solution 2/1 required changing the initial phase of the driving pressure. The quadratic model is unable to reproduce the subharmonic peak near 
$\nu /\nu _0 = 2$.
The third row of figure 2, panels (c) and (g), shows results for 
$R_0 = 100\ \mathrm {\mu }$m with 
$p_a/p_\infty = 0.5$. In the region of the main resonance, the biquadratic model appears to behave slightly better than the simpler quadratic one. However, neither model is able to follow the peak all the way to the amplitude of the Navier–Stokes simulation. Unlike the biquadratic model, the quadratic model is able to reproduce the subharmonic peak, but it proves distinctly less accurate below 
$\nu /\nu _0 \simeq$ 0.5. The biquadratic model does, however, reproduce the subharmonic peak if the driving is increased from 
$p_a/p_\infty = 0.5$ to 0.7 as shown later in figure 6.
The fourth row of figure 2, panels (d) and (h), shows results for 
$R_0 = 500\ \mathrm {\mu }$m with 
$p_a/p_\infty = 0.5$. Here we see a large gap in the region of the main resonance where even the Navier–Stokes solution fails to reach a steady state. At sufficiently large acoustic pressures, bubble oscillations can become chaotic (see e.g. Lauterborn & Kurz Reference Lauterborn and Kurz2010), a tendency that increases with the bubble radius and which is possibly the cause of these gaps. The sensitive dependence of chaotic oscillations on model details rules out a meaningful comparison of the Navier–Stokes solution with the approximate models. Aside from the main resonance, the biquadratic model performs reasonably well and proves clearly superior to the quadratic one.
While the graphs of 
$R_M^*$ versus frequency help to convey an overall picture of the model performance, they do not give any information on the precise waveform 
$R(t)$ of the oscillation versus time. Figure 3 shows a number of examples of this type for 
$R_0 = 10$, 100 and 500 
$\mathrm {\mu }$m in the first, second and third rows, respectively, all for the same pressure amplitudes as in figure 2. The solid lines are the Navier–Stokes solutions and the dashed lines the biquadratic model predictions. The examples shown have been chosen to be close to the peaks of the response curves in figure 2 and therefore illustrate, in some sense, the model performance in relatively difficult cases. Hardly any difference between the two sets of results can be seen for the 
$10\ \mathrm {\mu }$m bubble and the first three examples of the 100 
$\mathrm {\mu }$m bubble. The last example for this case is for a frequency in the subharmonic region that the approximate model cannot reproduce, and the difference with the Navier–Stokes solution is quite large. Differences between the Navier–Stokes and approximate models are more evident for the largest bubble (last row) but are modest.
Several examples of steady radius–time curves comparing the Navier–Stokes solution (solid lines) with the biquadratic model prediction (dashed lines). The four panels in the first row are for 
$R_0 = 10\ \mathrm {\mu }$m, 
$p_a/p_\infty = 0.8$ with 
$\nu /\nu _0 = 0.29$, 0.41, 0.72 and 1.96, respectively; in the second row for 
$R_0 = 100\ \mathrm {\mu }$m, 
$p_a/p_\infty = 0.5$ with 
$\nu /\nu _0 = 0.312$, 0.435, 0.90 and 1.95, respectively; and in the last for 
$R_0 = 500\ \mathrm {\mu }$m, 
$p_a/p_\infty = 0.5$ with 
$\nu /\nu _0 = 0.305$, 0.47, 1.2 and 0.237, respectively.
Going deeper into the model performance, we show in figure 4 the gas temperature distribution corresponding to the instants marked by circles in some of the panels of the last two rows of figure 3. We focus on the relatively more demanding cases of 
$R_0 = 100$ (first two rows) and 
$500\ \mathrm {\mu }$m, because, for smaller bubbles, very little difference is found. The solid lines (black) are the Navier–Stokes result, the dashed lines (red) the biquadratic model and the dash-dotted lines (blue) the adiabatic temperature at the same value of the radius (not necessarily reached at the same instant of time). The first row of figure 4 is for the first panel in the second row of figure 3 and it shows a very good agreement between the actual and biquadratic temperature distributions at the points of maximum radius, while the steepness of the distribution near the wall at the points of minimum radius, which approaches a boundary layer structure, is not captured as accurately. The effect of this error on the radius, however, is small, as can be seen from figure 3. The distributions in the second row of figure 4, corresponding to the third panel in the second row of figure 3, exhibit the same problem during the compression phase and can describe only in a general sense the turning of the temperature from hot to cold in the middle of the bubble expansion. The picture is not very different in the third and fourth rows of figure 4 for 
$R_0 = 500\ \mathrm {\mu }$m. Unlike the biquadratic distribution, the quadratic one cannot have a non-monotonic behaviour and the differences are larger. In most cases, the results given by the adiabatic model are significantly different.
Example of the gas temperature distribution at the instants marked by circles in figure 3. The first two rows are for 
$R_0 = 100\ \mathrm {\mu }$m with 
$\nu /\nu _0 = 0.312$ and 0.90; and the last two rows are for 
$R_0 = 500\ \mathrm {\mu }$m with 
$\nu /\nu _0 = 0.305$ and 1.2.
We show the model's performance in a demanding case characteristic of sonoluminescence by focusing on an example inspired by figure 4 of Brenner, Hilgenfeldt & Lohse (Reference Brenner, Hilgenfeldt and Lohse2002), i.e. the oscillations of a 
$4.5\ \mathrm {\mu }$m radius air bubble driven at 26.5 kHz by a pressure amplitude of 1.2 atm. For this case, liquid compressibility effects are important and we used the compressible version of the Rayleigh–Plesset equation shown in equation (20) of Brenner et al. (Reference Brenner, Hilgenfeldt and Lohse2002). This is the only change to our model made for the simulations of this case. The results for the radius–time behaviour are shown in figure 5, where the lines showing the Navier–Stokes solution and the biquadratic model prediction are hardly distinguishable from each other.
Comparison of the radius–time Navier–Stokes and biquadratic approximate solutions for a 
$4.5\ \mathrm {\mu }$m radius air bubble driven at 26.5 kHz by a pressure amplitude of 1.2 atm.
The conclusion from these and other results we have obtained is that the model works very well for smaller bubbles, for which the gas temperature distribution can be captured adequately by the biquadratic ansatz (3.3), while, in some cases, its performance may not be as good near resonances as the radius and driving pressure increase. In particular, the present approximation can be used to advantage in the modelling of bubbles in medical applications of ultrasound, for which typical bubble sizes are of the order of micrometres (see e.g. Dollet, Marmottant & Garbin Reference Dollet, Marmottant and Garbin2019; Helfield Reference Helfield2019; Versluis et al. Reference Versluis, Stride, Lajoinie, Dollet and Segers2020). As for the subharmonic response, which is important in these applications, we show in figure 6 that it can be reproduced by the biquadratic model for a sufficient driving pressure. Differences with the Navier–Stokes solution are to be imputed to the imperfect evaluation of damping on which the bubble behaviour in this frequency range is strongly dependent.
Detail of the response curve of a 
$100\ \mathrm {\mu }$m radius bubble driven at 
$p_a/p_\infty = 0.7$ in the subharmonic region, to be compared with the same bubble driven at 
$p_a/p_\infty = 0.5$ in figure 2(c).
5. Thermal effects
We now turn to a brief consideration of some further thermal aspects of the problem at hand. The three panels in figure 7 refer to a 
$100\ \mathrm {\mu }$m radius bubble driven by a pressure amplitude 
$p_a/p_\infty = 0.5$ at a frequency 
$\nu /\nu _0 = 0.9$. The radius–time behaviour for this case is shown in the third panel of figure 3 and some temperature distributions in the second row of figure 4; it is seen in figure 2 that, for this case, 
$R_M^* \simeq$ 0.9. Figure 7(
$a$) compares the centre temperature in the course of a steady oscillation as given by the Navier–Stokes solution (solid black line), the approximate biquadratic model (dashed red line) and the adiabatic model (dash-dotted blue line). The horizontal dotted line is the undisturbed liquid temperature and corresponds therefore to the isothermal model. The approximate model reproduces well the temperature at the centre of the bubble. The adiabatic model differs in several significant respects. In the first place, the maximum is shifted in time, which is a consequence of the insufficient damping of this model. Secondly, as a consequence of the same deficiency, the width of the peak is much narrower. This would cause larger radial accelerations, which may lead to the prediction of the fragmentation of the bubble by the well-known instabilities affecting the radial oscillations (see e.g. Plesset & Prosperetti Reference Plesset and Prosperetti1977). Finally, the gas is predicted to become much colder, with consequences that will be seen shortly.
Some results for 
$R_0 = 100\ \mathrm {\mu }$m, 
$p_a/p_\infty = 0.5$ and 
$\nu /\nu _0 = 0.9$ during one steady cycle: solid line (black), Navier–Stokes; dashed line (red), biquadratic model; dash-dotted line (blue), adiabatic; and dotted (purple), isothermal. (a) Centre temperature; (b) gas pressure versus bubble volume; and (c) normalized heat flow rate 
$Q^*=Q(t)/(8{\rm \pi} ^2 \nu _0 p_0 R_0^3)$.
Figure 7(
$b$) shows the bubble internal pressure versus volume during a cycle. The Navier–Stokes and approximate results are difficult to distinguish and nearly overlap, enclosing a thin region the integral of which (when plotted on linear scales) is the energy dissipated over a cycle. The area of this loop appears small because the energy loss per cycle is small compared with the mechanical energies at play (see figure 8 below), but nevertheless this energy loss is much greater than that due to viscosity or (at this pressure) acoustic radiation. The dashed blue and dotted purple lines are the adiabatic and isothermal results, respectively. In view of the functional dependence of 
$p$ upon 
$V$ that these models presuppose, these lines can only enclose a zero area. It is seen that the adiabatic model predicts a much smaller minimum volume and a larger maximum volume than the Navier–Stokes solution. Owing to the lack of heat exchange with the liquid of this model, the temperature and, therefore, the pressure are too low during the expansion. Accordingly, upon compression, the pressure starts increasing too late and is therefore unable to slow down the inward motion sufficiently fast. The very large pressure ultimately reached then causes a larger expansion of the bubble beyond the correct value. For obvious reasons, the isothermal minimum volume is even smaller than the adiabatic one. Since, as (1.2) shows, the viscous dissipation is inversely proportional to 
$R$, the very small radii predicted by the adiabatic and isothermal models enhance the effect of viscosity much above its proper level.
Normalized rate of work 
$\dot {W}^*=\dot {W}/(8{\rm \pi} ^2\nu _0 p_0R_0^3)$ (outer loops) and heat flow rate 
$Q^*$ (inner loops) versus bubble radius in the course of a steady oscillation. Solid lines (black), Navier–Stokes; dashed lines (red), biquadratic model. (a) 
$R_0 = 100\ \mathrm {\mu }$m, 
$p_a/p_\infty = 0.5, \nu /\nu _0 = 0.9$; (b) 
$R_0 = 10\ \mathrm {\mu }$m, 
$p_a/p_\infty = 0.8, \nu /\nu _0 = 1.1$; (c) 
$R_0 = 1\ \mathrm {\mu }$m, 
$p_a/p_\infty = 1.6, \nu /\nu _0 = 0.92$.
It is interesting to note that, in the latter part of the expansion phase, the 
$p$–
$V$ behaviour is nearly isothermal, while it is closer to adiabatic in the compression phase due to the slowness of the former and to the rapidity of the latter. A feature worthy of note is that the loop is displaced above both the isothermal and adiabatic lines. The reason is rectified heat transfer into the bubble, which causes the gas mean temperature to increase until the extra heat loss compensates the heat rectification.
Figure 7(c) shows, as a function of time during a cycle, the instantaneous heat transported out of the bubble given by 
$Q(t)=4{\rm \pi} R^2 [k \partial _rT]_{r=R}$, with the solid and dashed lines representing the Navier–Stokes and model predictions. The agreement is far from perfect, although the model captures the essence of the behaviour.
Another way to represent the thermo-mechanical behaviour of the bubble is by showing the rate of work 
$\dot {W}=-p\,\textrm {d} V/\textrm {d} t$ and heat flow 
$Q(t)$ versus radius as in figure 8 for, from left to right, 
$R_0 = 100\ \mathrm {\mu }$m, 
$p_a/p_\infty = 0.5, \nu /\nu _0 = 0.9$ (same as in the previous figure), 
$R_0 = 10\ \mathrm {\mu }$m, 
$p_a/p_\infty = 0.8, \nu /\nu _0 = 1.1$, and 
$R_0 = 1\ \mathrm {\mu }$m, 
$p_a/p_\infty = 1.6, \nu /\nu _0 = 0.92$. The solid lines are the Navier–Stokes solution and the dashed lines the biquadratic model predictions. The mechanical work 
$\dot {W}$ is seen to be predicted better than the thermal energy loss, although the latter is much closer to the correct result for the smaller 10 and 
$1\ \mathrm {\mu }$m radius bubbles. These results are another confirmation of the fact that the predictions of the approximate models improve significantly as the bubble radius decreases, as the wall heat flux tends to be less extreme for smaller bubbles.
6. Conclusions
We have presented two approximate models for the thermal interaction of a gas bubble with the host liquid. The comparison with a first-principles Navier–Stokes simulation of the thermo-mechanical process has shown the models to be accurate and therefore useful over a broad range of parameters. In other parameter ranges, the incorrect estimation of the energy loss to the liquid makes the models less accurate, but they still offer an alternative to the polytropic model, which neglects thermal losses entirely. Thus, in spite of their limitations, these models afford at a modest computational cost a more faithful representation of the physics of the bubble–liquid thermal interaction than currently possible without recourse to partial differential equations.
The correct use of the models requires an understanding of their limitations, which, for both, and for the quadratic model more than for the biquadratic one, lie in an inability to reproduce the temperature boundary layer that develops with larger bubbles, larger-amplitude oscillations and faster processes. This deficiency prevents the correct calculation of the heat exchange with the liquid. If the net heat loss is too small, the initial transient lasts too long to allow the bubble to reach the steady-state regime at the correct rate. This is what may prevent the approximate models from reproducing the large-amplitude results near the fundamental resonance in figures 2(c), (d), (g) and (h). Secondly, the point of vertical tangency of figure 1 is dependent on damping and, therefore, the approximate models may be inaccurate near the points in which a jump to the higher amplitude would occur. Thirdly, some features of the response, such as the subharmonic, are strongly dependent on the correct amount of damping and, indeed, the approximate models perform unevenly in some of these cases. The model of Toegel et al. (Reference Toegel, Hilgenfeldt and Lohse2002) (see also Stricker et al. Reference Stricker, Prosperetti and Lohse2011) incorporates a boundary layer structure and may be preferable in the parameter ranges where the models presented here have deficiencies.
Acknowledgements
This research was made possible by a grant from the Gulf of Mexico Research Initiative. Data are publicly available through the Gulf of Mexico Research Initiative Information & Data Cooperative (GRIIDC) at https://data.gulfresearchinitiative.org (doi: 10.7266/M4DHDRBG). The numerical computations were carried out on the Sabine cluster of the University of Houston Research Computing Data Core.
Declaration of interests
The authors report no conflict of interest.
