The most popular technique for estimating the gas bubble size distribution (BSD) in liquids is through the inversion of measured attenuation and/or sound speed of a travelling wave. The model inherent in such inversions never exactly matches the conditions of the measurement, and the size of the resulting error (which could well be small in quasi-free field conditions) cannot be quantified if only a free field code exists. Users may be unaware of errors because, with sufficient regularization, such inversions can always be made to produce an answer, the accuracy of which is unknown unless independent (e.g. optical) measurements are made. This study was commissioned to assess the size of this error for the mercury-filled steel pipelines of the target test facility (TTF) of the spallation neutron source at Oak Ridge National Laboratory, TN, USA. Large errors in estimating the BSD (greater than 1000% overcounts/undercounts) are predicted. A new inversion technique appropriate for pipelines such as TTF gives good BSD estimations if the frequency range is sufficiently broad. However, it also shows that implementation of the planned reduction in frequency bandwidth for the TTF bubble sensor would make even this inversion insufficient to obtain an accurate BSD in TTF.
The purpose of this paper is to devise a method of estimating the number and sizes of gas bubbles in liquid from the measured attenuation and/or phase speeds of a travelling acoustic wave that is suitable for the mercury-filled steel pipes of the target test facility (TTF) of the $1.4 billion spallation neutron source (SNS) at Oak Ridge National Laboratory (ORNL), TN, USA. A further aim is to use that method to test whether a sensor based on such principles is practical when it emits the frequency range originally designed, and when it emits the narrower but more affordable frequency range enforced by the 2008 global financial crash.
In an inversion (the estimation of a parameter from measurements where there is not a simple one-to-one mapping and a model is required), a basic ideal requirement is that the assumptions of the model used in the inversion match the conditions under which the measurements are made. Strictly speaking, this has never happened for the most widespread technique of measuring bubble populations acoustically, when the bubble size distribution (BSD) is estimated from the measured attenuation and/or sound speed of a travelling wave (Commander & McDonald 1991; Duraiswami et al. 1998; Leighton et al. 2004). This is because the model used assumes acoustic plane wave propagation in free field (PWFF) conditions (Commander & Prosperetti 1989).
Some other techniques do not require the measurement of attenuation and sound speed of travelling waves in bubbly liquid. Important examples are the increasingly popular resonator technique (Farmer et al. 2005) and the impedance tube method of Wilson et al. (2005). Both however could be affected by coupling between the liquid and the container walls, and indeed Wilson et al. carefully designed their impedance tube to avoid this regime, stating: ‘This impedance tube was designed such that when filled with bubble-free water, both the radial component of the plane wave mode and the higher-order mode are negligible and can be ignored…. An inversion scheme that properly accounts for both components of both modes does not exist and for now these elastic waveguide effects have been ignored’ (Wilson et al. 2005; p. 1898). They later used the lossless model of Lafleur & Shields (1995) to predict the effect of coupling on the phase speed in bubbly water, but could not predict attenuation or invert for the bubble population. Boston University collaborators on ORNL's bubble detection programme attempted to use a mercury-filled, flow-through, stainless steel cylindrical resonator to measure sound speed. Although their system worked well with water (Ormonde et al. 2008; Roy et al. 2008), they encountered problems with imperfectly wetted steel/mercury boundaries, particularly when bubbles were introduced (R. A. Roy 2012, personal communication).
Given coupling was a concern to Wilson et al. when studying propagation in water-filled tube with thick steel walls, it is not surprise that it becomes a dominant feature of propagation in the mercury-filled steel pipes of the TTF at ORNL. The current paper provides an inversion method which included this effect for measurement of attenuation and speed of travelling waves.
The mismatch between model and scenario has not prevented the PWFF inversion being used in the past, since with sufficient regularization it can always be made to generate an answer, and rarely is an independent measurement made to check the correctness of that answer. In some scenarios (e.g. with low bubble concentrations in deep oceans), the departure from PWFF conditions may produce insignificant errors (Vagle & Farmer 1998). However, the petrochemical, cement, food production, ceramic, power, mining and pharmaceutical industries, etc., require methods for measuring bubble populations in liquids in pipelines (Campbell & Mougeot 1999; Iskandrani & Kojasoy 2001; Punurai et al. 2006; Yim & Leighton 2010). Here, the error produced in the mismatch between the environment and the PWFF assumption used in the inversion has never been quantified, let alone corrected for. This paper provides these two advances. In §2, simulated data are used, because for this the input bubble population is known and the errors in the inversion are clearly quantified. It is shown that the PWFF inversion produces large errors in the estimation of the BSD (greater than 1000% overcounts/undercounts). A new method of inversion is provided, which takes into account the fact that the measurements were not made under PWFF conditions but instead occurred in a pipe where coupling occurs between the pipe walls and the fluid inside it. This inversion is validated against simulated data from pipes and is shown to give similar levels of accuracy to those obtained when the PWFF inversion is applied to data taken in PWFF conditions. Section 3 describes experimental tests on water-filled poly methylmethacrylate (PMMA) pipes, which had previously been shown to mimic the acoustic coupling that would occur in the mercury-filled steel pipes of ORNL's SNS TTF (Baik et al. 2010; Jiang et al. 2011). Although these data show that the new inversion technique would provide accurate BSDs if the frequency range of the original sensor design were to be used, budget cuts forced by the 2008 global financial crisis dramatically reduced the frequency range that could be afforded. Section 4 demonstrates that it is not worth ORNL commissioning a sensor with this reduced frequency range, because although it is affordable, the loss of accuracy is unacceptable.
Measurement of the bubble population is however important for ORNL, since helium bubble injectors are planned to be fitted to the SNS in 2013 to absorb the shocks generated when the proton beam impacts the mercury. Such absorption will mitigate cavitation damage to the casing surrounding the mercury (Chitnis et al. 2010; Manzi et al. 2010; Leighton et al. 2011).
Assessment of the effect of applying PWFF inversions to pipe data, and assessment of the investment which needs to be made to cover enough frequencies to produce satisfactory results, are vital steps if industry is to make use of acoustic inversions for BSD in pipes. This approach should replace the current faith that the technique can produce an answer (as it always can be made to do) without questioning the accuracy of that answer.
The theory underlying the PWFF inversion technique assumes that plane waves propagate, with a single sound speed and known frequency-dependent attenuation in bubble-free conditions. When bubbles are added in PWFF conditions, at any given frequency there exists a single sound speed (which is frequency dependent), and the attenuation changes. These changes are attributed wholly to the presence of bubbles. When modal propagation occurs in pipes where the liquid is coupled to the walls, even in bubble-free conditions each mode propagates with its own frequency-dependent longitudinal phase velocity and attenuation (Baik et al. 2010). The differences between these phase speeds and attenuations, and those assumed for infinite volumes of bubble-free liquid, would be attributed by the PWFF inversion as being caused by bubble presence. In such circumstances, if an industry cannot cancel out this error by generating bubble-free conditions and instead relies on a predicted bubble-free attenuation, then use of a PWFF inversion would give a finite bubble count even when no bubbles are present. The errors in the inversion when bubbles are added have not previously been investigated. It should be noted that there are further sources of mismatch that are not included in the formulation of this paper, but which have been considered elsewhere (§4a).
The methods developed in this paper have also been used to provide demonstrations for TV, public shows and student demonstrations of the ability of bubbles in liquids, and aerosol drops in air, to absorb sound (Leighton et al. 2011, 2012).
(a) The forward problem
For pipes with boundary conditions resembling those of the mercury-filled steel pipes of ORNL or the water-filled PMMA pipes of §3, Baik et al. (2010) showed how the complex axial wavenumber, q′0m of the axisymmetric coupled modes that propagate in the axial direction is obtained by their eqn (6) with the matrix elements of their appendix B (the ‘0’ in the subscript ‘0m’ refers to the mode being axisymmetric and ‘m’ refers to the mode index). Eqn (6) of that paper gives the infinite number of the axisymmetric modes (denoted as ETm where m is mode index related to the radial motion) and the real and the imaginary parts of the complex solution, q′0m, which give the phase speed and the attenuation of the modes respectively in a pipe containing a bubble-free liquid.
The introduction of bubbles into the previously bubble-free liquid in the pipe changes the elastic properties of the liquid. The complex wavenumber k′1 in the liquid is far more influenced by the bubble-induced changes in compressibility than it is by the changes in the spatially averaged density of the bubbly mixture (at least up to gas volume void fractions of 10%). Changes in the complex wavenumber k′1 affect both phase speed and attenuation, and bubbles moreover scatter the sound and change the stress elements and displacement vectors in the liquid contacting the inner pipe wall.
The sound velocity of a bubbly liquid in a thin elastic tube has been calculated using a one-dimensional approach and the balance between equations of mass, momentum and energy (Rath 1981; Schubert 1991). The result predicts subsonic sound speeds as a function of the gas volume void fraction, but it is not sufficient for the current study because acoustic coupling between the elastic solid and the bubbly liquid will in practice also give an infinite number of supersonic modes. To encompass such modes in this study, the characteristic equation for the axisymmetric modes in a tube filled with bubbly liquid is implemented by combining the description of acoustic modes in a pipe containing bubble-free liquid (Baik et al. 2009, 2010) with an analysis which predicts the bubble-induced changes in phase speed and attenuation that occur in an infinite bubbly liquid (Commander & Prosperetti 1989). Wilson et al. (2005) applied the formulation of Commander & Prosperetti (1989) in order to measure the phase speed and attenuation that occurred in an impedance tube made of heavy stainless steel, where the phase speed of the lowest axisymmetric mode at low frequencies does not differ significantly from the intrinsic sound speed of the liquid within the tube. In their approach, the phase speed of the dominant axisymmetric modes in the bubble-free waveguide is used as an input parameter. In contrast, the current study uses the BSD as a parameter to obtain the dominant axisymmetric modes in the bubbly waveguide.
For a given homogeneous size distribution of bubbles in an infinite liquid undergoing steady-state pulsations in the linear long-wavelength limit (Clarke & Leighton 2000; Leighton et al. 2004), the complex phase speed of sound in bubbly liquid, Cb, is given, with the −iωt convention, by 2.1where C1 is the intrinsic sound speed in bubble-free liquid (Commander & Prosperetti 1989), R0 is the equilibrium bubble radius, ω0 is the resonance frequency and b is the total damping parameter (Ainslie & Leighton 2011). The presence of bubbles changes the acoustic impedance of liquid, reducing it at low frequencies in infinite media (Mallock 1910; Leighton & Robb 2008), but tending to increase it at frequencies just above the bubble resonance (Commander & Prosperetti 1989). This change is described by the transformation of the real wavenumber in bubble-free liquid, k1, into the complex wavenumber in bubbly liquid, k′1. The BSD is described by nb(R0), where nb(R0)dR0 is the number of bubbles per unit volume having radii between R0 and R0+dR0 (where conventionally dR0=1 μm when such populations are represented graphically). The gas volume void fraction, Γ, is given by . The total damping parameter, b, includes bth, bvis and brad, corresponding to thermal, viscous and radiation damping, respectively, and explicit forms of ω0 and b are given by Commander & Prosperetti (1989) (those authors, who unlike here write with the+jωt convention, have identified and thoroughly publicized a typographical error in their expression for bth; Ainslie & Leighton 2011). A key assumption is that the complex wavenumber k′1 in equation (2.1) refers to the wavenumber in the infinite body of the bubble-free liquid. This must be adapted for the situation when the bubbly liquid is contained within a tube, where the bubbles change the pressure field within the pipe and, accordingly, the complex wavenumber k′1 also changes. It is assumed that the bubbly liquid axisymmetrically exerts on the tube wall an average pressure which is established by equation (2.1). Therefore, when a bubbly liquid with a known BSD fills an elastic pipe, the axisymmetric modes within the pipe are calculated by substituting the expression of k′1 from equation (2.1) into eqn (6) (with the matrix elements of appendix B) of Baik et al. (2010).
This is done in figure 1. The solid curves predict the phase speed (identical solid curves in (a) and (c)) and attenuations (identical solid curves in (b) and (d)) of three modes when the PMMA pipe is filled with bubbly water. The BSD of the bubbly water is based on the smoothed version of the optically measured BSD (the solid black μCORT curve of figure 3), where equation (2.1) in this paper and eqn (6) of Baik et al. (2010) are used to predict all the solid curves in (a–d). These are plotted as a function of wavenumber-radius product, x=k1a, where a is the inner radius of the tube. To show the effect made by the addition of bubbles, the dashed lines in figure 1a,b show the modal phase speeds and attenuations when the same pipe is filled with bubble-free water. To show the effect made by the pipe walls, the dashed lines in figure 1c,d show the modal phase speeds and attenuations for PWFF propagation in an infinite (unconfined) volume of bubbly water with the same BSD (as calculated by equation (2.1)). This clearly indicates the degree of error which would occur if the measured sound speed in a pipe were inverted using PWFF theory to infer BSD, the method currently used in such inversions (Hansen et al. 2004). Looking at figure 1a for example, an accurate inversion should be based on this change in the modal speeds caused by the introduction of bubbles. It should not be based (as the PWFF inversions would do) on the departure of any curve from the horizontal line given by C0m/C1=1, an approach made worse when the ‘time of first arrival’ is used to calculate sound speed (without recognizing modal propagation) as this selects data for the faster modes.
If the only measurements available were the dispersion relation of the modes with (solid curves in (a) and (b)) and without (dashed curves in (a) and (b)) bubbles present, could the BSD be inferred uniquely? The two sets of dispersion curves are sufficiently similar to allow identification of features, but the question is whether these features change sufficiently to act as a good indicator of BSD directly. For example, the break points (Baik et al. 2010) where two modes come into closest approach (barring the tendency to converge as ) occur in bubble-free conditions at x∼2.8 for ET1 and ET2, and x∼6.4 for ET2 and ET3. When this BSD is introduced, these break points shift to slightly higher frequencies. Another potential indicator is that the phase speeds of the supersonic modes in the bubbly pipeline are higher than those in bubble-free conditions over most of the frequency range shown.
Figure 1b shows that throughout most of this frequency range in the pipe, for this BSD the attenuation of a given mode in bubbly water is larger than that in the bubble-free case. This might be expected from the high absorption per metre seen in infinite volumes of bubbly liquid (shown for convenience using the dashed curves in figure 1d). It is interesting to note that for ET1 the attenuation rate in the pipe is less (e.g. in the frequency range of x<4) than for plane waves in an infinite body of water for the same BSD as shown in figure 1d. The reason for this is that the longitudinal and the shear ultrasonic absorptions in the PMMA are less than the absorption in free field bulk bubbly water in that frequency range. As a result, the acoustic coupling between bubbly water and the PMMA results in the intermediate attenuation of the ET1 mode between the absorptions of the PMMA and bubbly water.
In summary, it is important to note that the significant discrepancy between the phase speeds and attenuation for plane waves in an infinite body of liquid (the dashed curves in figure 1c,d) and the values predicted for the same BSD contained within the pipe (solid lines) means that inversions that apply the standard methods (Commander & McDonald 1991; Duraiswami et al. 1998) to bubbly water contained in vessels will generate errors in the inferred BSD.
(b) The inverse problem
It is now routine to compare the attenuation measured across a wide range of frequencies between two hydrophones before and after the injection of bubbles to infer the bubble population through inversion in conditions of low void fraction (Commander & McDonald 1991). Although in principle this can also be done for phase speed (Duraiswami et al. 1998), in practice, this is not so useful, since small changes in phase calibration of the hydrophones, or displacements or inaccuracies in hydrophone position, can give large artefacts, and at the very least good practice would require that the measurement be repeated with the hydrophone positions swapped to an accuracy very much less than the smallest wavelength used.
Even if used in actual PWFF conditions, there are complications with such inversions: for example, the algorithms assume no bubbles exist that are resonant outside of the bandwidth of insonification, so that attenuation caused by any bubbles that occur in these two forbidden size regions must be accounted for by erroneously enhancing the bubble count within the allowed bandwidth. An iterative approach to estimating the BSD can be undertaken to mitigate this effect, but this is rarely done (Caruthers et al. 1999). Most importantly, with sufficient regularization almost any measured attenuation could be inverted to give an estimated BSD, but there is no in-built validation that this estimate of the BSD is accurate. This can be illustrated by the way in which changes to the regularization can yield different estimates of the BSD for the same input data. Hence provision of an answer from this method does not indicate that it is the correct answer.
Assume that the liquid is effectively infinite and homogeneous, that the void fraction is ≪∼0.0001 unless multiple bubble–bubble interactions are taken into account (Kargl 2002) and that all bubbles are spherical and pulsate linearly in the steady state (Clarke & Leighton 2000). If the perturbations in the radius of the oscillating bubble, ΔR, are small (i.e. ΔR≪R0), then from equation (2.1) the real and imaginary parts of the squared complex wavenumber, , can be respectively approximated as follows 2.2where uψ=Re[k′1]/k1 and υψ=Im[k′1]/k1 are the ψth elements of the normalized real and imaginary parts of the complex wavenumber k′1. The subscripts ξ and ψ are necessary in the statement of ω0ψξ because of the dependence of the resonance frequency in the matrix element on the driving frequency through the frequency dependence of the polytropic index (Ainslie & Leighton 2011). Because M is the total number of discretized elements in the implementation of equation (2.2), the sampling interval of bubble radius ΔR0 is related to M by ΔR0=(R2−R1)/(M−1) where R1 and R2 are the minimum and maximum equilibrium radii contained in the bubble distribution. Equation (2.2) shows how uψ and υψ are related to each other for a given frequency element ωψ. Stacking equation (2.2) along a frequency range of [ω1,ω2] (where ω1 and ω2 are minimum and maximum frequencies in the measurement) recasts equation (2.2) into the following matrix form: 2.3the matrix 1 is a N×1 matrix representing unity throughout a given frequency element ωψ, and Nb is a M×1 matrix of which each element is the bubble size spectrum nb(R0ξ) at a given equilibrium radius of bubble, R0ψ. The matrix Z has a dimension of N×M and each element Zψξ can be represented as 2.4where Yψ=υψ/uψ is the ratio of the imaginary part to the real part of the complex wavenumber k′1. The bubble resonance frequency, ω0, and total damping parameter, b, are also functions of ωψ and R0ξ (Ainslie & Leighton 2011). Equation (2.4) requires knowledge of both the phase speed change and the attenuation caused by bubbles in order to estimate the BSD, Nb. Such knowledge is provided by the measurement of Yψ=υψ/uψ in the formulation. The conventional acoustic inversion problem (which only requires measurement of the attenuation) is reconstructed when only the term associated with 1/Yψ inside the square bracket of equation (2.4) is considered, and the condition of Re[k′1]=k1 is imposed. Estimation of the measured BSD, Nb, from the attenuation (α) that was measured across a sufficiently wide frequency range (i.e. one that encompasses the pulsation resonances of all bubbles present) is a problem of linear algebra that is an ill-posed inversion requiring regularization (Tikhonov & Arsenin 1977; Commander & McDonald 1991; Hansen 1998; Leighton et al. 2004). The current study uses Tikhonov regularization (Tikhonov & Arsenin 1977; Hansen 1998), which requires the user to define a regularization parameter close to the optimal value. In order to find the optimal value, the current study adopted the L-curve method which finds the vertex of the curve that plots the Euclidian norm of the regularized solution against the norm of the error (Leighton et al. 2004).
As stated in §2a, PWFF theory assumes that there exists a single frequency-dependent attenuation, and a single sound speed that is frequency-independent in the bubble-free condition and frequency-dependent when bubbles are added. However, in pipelines both sound speed and attenuation are multi-valued at any given frequency and are functions of frequency in both bubbly and bubble-free conditions. The two sets of observations cannot match over a broad frequency range, as was shown by the mismatch between the dashed and solid lines in figure 1c,d. Application of equation (2.3) to invert attenuation to estimate BSD in pipelines would be erroneous. There are further sources of mismatch that are not considered by this paper (§4a).
In order to invert the complex wavenumber measured in a bubbly pipeline to estimate the BSD, this paper develops a technique based on estimating the complex wavenumber that would occur in the infinite body of bubbly liquid that had the same BSD as was present inside the pipeline. This is achieved by applying a Taylor expansion to the formulation of the forward problem described in §2a. For a given complex wavenumber (k′1) in the liquid contained inside the tube, the complex wavenumber of the propagating mode inside the tube, q′0m (the real part of which is used to calculate C0m) satisfies the following equation: 2.5where F appears in Baik et al. (2010) as the left-hand side of their eqn (6) with the matrix elements in their appendix B. Let the complex wavenumbers in the infinite body of bubble-free liquid and bubbly liquid be k′1f and k′1b, respectively (where the infinite bubbly liquid has the same BSD as is occurring in the pipe). Let the complex wavenumbers of the axisymmetric modes in the bubble-free liquid and bubbly liquid inside the pipe be q′0mf and q′0mb, respectively. Suppose that the phase speed and the attenuation are measured in the pipe when it is filled with the bubbly liquid. This measurement gives q′0mb. The complex wavenumbers, k′1f and q′0mf, can be measured, although it is safest to check that nominally bubble-free measurements agree with the truly bubble-free predictions of Baik et al. (2010) in cases where, as in TTF, the convoluted sealed pipework has a history of previous injection of barely soluble gas. The only unknown variable is the complex wavenumber, k′1b, which is the input parameter in the acoustic inversion problem to estimate the BSD. This is obtained from the following step. The variables, k′1f, k′1b, q′0mf and q′0mb should satisfy equation (2.5), that is: 2.6and 2.7Provided that the introduction of bubbles inside the pipe does not change phase speed and the attenuation of the modes much from those in bubble-free liquid such that |q′0mb/q′0mf−1|≪1 (figure 1), then equation (2.7) can be expanded to the first order by Taylor expansion: 2.8The current study adopted the model by Commander & Prosperetti (1989) to describe the complex wavenumber in a bubbly liquid that was developed from the linearized Keller equation. We assume that the small deviation between axial wavenumbers in bubble-free liquid and in bubbly liquids is linear, allowing use of a Taylor expansion to map from one to the other. Although the choice of a linear method imposes technique-sensitive limitations, it has the advantage over other options (such as use of Bayesian statistics) in the situation where measurements are only possible over certain limited frequency regimes. In such circumstances, the probability distribution (defined in terms of frequency or number of modes) does not change smoothly, generating discontinuities in the probability distribution.
From equations (2.6)–(2.8), the changes in the complex wavenumber in liquid are 2.9Therefore, k′1b is estimated to be k′1f+Δk′1, and this is used as an input parameter of the acoustic inversion described in equation (2.3). The bubble-free and bubbly conditions inside the tube are accounted for by Δq′0m which contains the change of phase speed and additional damping.
Numerical error is always a problem when equation (2.9) is applied to an inverse problem to estimate the BSD contained in the pipe. Such errors contain values of k′1b which are either too small or too large. In the worst case, the value of k′1b estimated through equation (2.9) generates negative values, which will be shown in figure 4a. These numerical errors come from steep changes of derivatives of the characteristic equation F(k′1,q′0m) with respect to change of q′0m or k′1. Sometimes, it changes phase, and therefore negative elements of the wavenumber can be obtained within the linear variances of q′0m and k′1. Of course such errors can be mitigated by expanding the Taylor series to higher order. However, this makes the estimation of k′1b more complicated, and even with higher-ordered expansions, such numerical errors still exist whenever large changes between bubbly and non-bubbly conditions are observed. Therefore, in the actual implementation of equation (2.9), limitations of the range of Δk′1 are imposed so that if the resulting estimate of k′1b were to be too large or too small, it is abandoned and q′0mb is used in its place. It is apparent in figure 1 that the differences between the imaginary parts (corresponding to the attenuation) of k′1b (for an infinite volume of bubbly water) and q′0mb (for bubbly water in the pipe) are greater than those between the real parts of k′1b and q′0mb across most of the frequency range. This is due to the fact that bubbles change the attenuation more effectively than they do the phase speed of media. Thus, in the present study, the extent of possible values of k′1b was set to go from 0.5Re[q′0mb] to 2Re[q′0mb] in the real part, and from 0.1Im[q′0mb] to 1000Im[q′0mb] in the imaginary part, ranges which were estimated from the data in figure 1. Outside of these ranges, k′1b was simply substituted with q′0mb.
3. Experimental method
Experiments were carried out in a water-filled PMMA tube of 2 m length, 4.445 cm inner radius and 0.5 cm wall thickness which is supported vertically (Baik et al. 2010; Leighton et al. 2012). The driving frequency of the function generator, which does not necessarily equal the peak frequency of the energy in the water (Baik et al. 2010) was incremented from 15 to 35 kHz in 1 kHz steps, with 20 cycles per pulse. To obtain the sound speed and attenuation at each frequency, the two-dimensional Fourier Transform technique and Prony's method were adopted (Baik et al. 2010; Jiang et al. 2011).
Air was injected 15 cm above the tube base using a lumbar puncture needle (90 mm length and 0.6 mm outer diameter), the needle being vibrated using a mobile phone vibrator to prevent bubble coalescence at the needle tip (Leighton et al. 2011, 2012), changing the bubble radius range from 1–2.5 mm to approximately 30–1000 μm.
An optical method, called μCORT (Multi-Color Rise Time, a version of which also exists in monochrome—Leighton & White 2012) provided an independent (non-acoustical) estimate of the BSD by exposing each frame of the cloud of rising bubbles to a red flash, followed by a blue flash a short time later (0.025 s for figure 2; 0.042 s for the 0.0058% void fraction of figure 3). This removes the need to image the bubble wall accurately (and assume a given dimension in the third dimension) that are requirements of most photographic techniques for obtaining BSD, and therefore restricts them to the bubble size range that can be imaged accurately in the field of view. By eliminating that need, μCORT overcomes such limitations, expands the field of view and depth of field, and extends the bubble radius range. μCORT is immune to image distortion caused by the curvature of the pipe walls, and overcomes errors introduced by software that cannot recognize that what it might perceive to be a bubble perimeter in fact only extends to a highlight of the illumination on a fraction of the bubble wall. μCORT automatically calculates the rise speed of each of the thousands of bubbles in the field of view by correlating the locations of the neighbouring red and blue images in the frame to pair up automatically the images of the same bubble, and hence its displacement during the interval between the flashes. The potential errors come in then converting those rise speeds to bubble radii. While calibration curves exist which might allow conversion of rise speed into bubble size for air–water mixtures such as used here (Clift et al. 1978), the tendency of the bubbles to act collectively, and for the pipe geometry to affect rise speed (Krishna et al. 1999; Mukundakrishnan et al. 2007), must be considered. The rise times of bubbles in swarms in tubes of similar size to the one used here have been studied previously (Garnier et al. 2002; Simonnet et al. 2007). Given that Garnier et al. (2002), working with larger bubbles (up to 2.75 mm radius) at higher void fractions (1–40%) than used here (approx. 30–1000 μm and 0.005%, respectively) concluded that wall effects could be neglected in terms of rise speed, and that the results were applicable for void fractions of less than 1 per cent, that assumption is also applied here. However, Garnier et al. (2002) indicate that, for bubbles smaller than 2.75 mm radius, if an isolated bubble were to rise at a speed in a quiescent liquid, then its equivalent rise speed (Vrel) relative to the liquid flow when it is part of a cloud of void fraction Γ is . This correction was tested by Guet et al. (2004) for bubbles smaller than 3 mm radius and void fractions smaller than 20 per cent (scenarios which cover the current case). This correction, though small (approx. 4%) for the void fractions (less than 10−4) studied here, is applied in this paper.
The effect of containment in a pipe on bubble resonance is rarely calculated (Leighton et al. 1998a, 2002). Using the worst case formulations of Leighton (2011), neglect of this effect causes a systematic error which underestimates the bubble inertia associated with its pulsation resonance by 3 per cent for bubble radii of 30 μm, rising to approximately 100 per cent for bubble of radius 1mm. However, the lowest frequency applied to the population is here 15 kHz, corresponding to a resonance with a bubble of radius approximately 200 μm where the inertia is underestimated by 20 per cent. Though significant, an iteration step was not applied to undertake this correction because it was small compared with the error discussed in §4.
The grey solid curve in figure 3 shows the raw optical BSD (bubble population present during the acoustic tests) obtained by the μCORT method. In order to obtain smooth curves of calculations for bubbly media, raw optical BSD is digital filtered to produce the black solid curve, which is the BSD adopted throughout this paper. The meanings of other symbols (circles, triangles, squares and asterisks) are explained in the figure caption and §4a. Note that the average bubble radius is approximately 300 μm, and for such bubbles the resonance frequency (approx. 10 kHz) is lower than the lowest frequency the affordable acoustical source can generate, illustrating the practical limitations of acoustic inversions of this type. This BSD has a gas volume void fraction of Γ≈0.0058%.
(a) Validation tests using simulated acoustic data based on optical measurements
The optical (μCORT) estimation of the BSD that occurred in the pipe during the acoustical tests (the black solid line in figure 3) was used as input in the forward model of Commander & Prosperetti (1989) (equation (2.1) in this paper) to predict what phase speed and extra attenuation the bubbles would have generated for PWFF propagation in an infinite body of water (the dashed curves in figure 1c,d). When these dashed line data in figure 1c,d are used as input in an inversion which assumes an infinite body of bubbly liquid following the method of Commander & McDonald (1991), equation (2.3), that inversion of course recovers the BSD used as the original input (the open circles in figure 3, with characteristic discrepancies at the ends of the size range illustrating the error introduced by the inversion itself). Such agreement is of course meaningless to the circumstances that actually exist in the pipe.
However, the importance of figure 3, and the raison d’être of this paper, is in showing what would happen if sound speed and attenuation data taken in a pipe were to be used to estimate the BSD under the traditional assumption of PWFF conditions. The acoustic data are simulated (the ET1, ET2 and ET3 solid curves in figure 1, calculated using the theory of Baik et al. 2010) because that allows us to know the actual BSD data used as input (solid black line in figure 3), against which we can compare the BSD estimated using PWFF theory in the traditional manner. The inversion proceeds as follows. The predicted complex wavenumber of the ET1, ET2 and ET3 modes for bubbly conditions (the solid curves in figure 1) are directly used as simulated input data Yψ. These data are then inverted using equation (2.3) for the inversion (which assumes an infinite body of liquid). Since the in-pipe conditions under which the solid line ET1, ET2 and ET3 acoustic data of figure 1 where ‘measured’ are not matched by the PWFF assumptions of equation (2.3), the estimates of the BSD are expected to be wildly erroneous. This is indeed shown to be the case in figure 3 for ET1 (triangles), ET2 (squares) and ET3 (asterisks). This demonstrates that the traditional inversion technique, uncorrected for deviations from PWFF conditions, cannot be assumed to be accurate when such deviations occur. For example, the ET1 data (triangles) above approximately 200 μm radius in figure 3 record more than 1000 per cent undercounts, while the other symbols (squares and asterisks) record more than 1000 per cent overcounts (too large to be plotted within the window of the figure). Therefore, the standard PWFF inversion can produce large errors in the estimation of the actual BSD inside the pipeline.
There is, however, a solution. As explained by equation (2.9), the complex wavenumbers that would have occurred were the same BSD present in PWFF conditions can be inferred from figure 1's simulated data for the phase speed and the attenuation of the modes in a pipe containing that same BSD. This can be used to ensure that the inversion is undertaken on data which complies with the assumptions inherent in that inversion. Since PWFF conditions never exist in practice in bubbly liquids, strictly speaking no such inversion for BSD from travelling wave attenuations and sound speed has ever done this, although clearly under some circumstances the approximation is sufficiently good (Wilson et al. 2005). Note however that the method here does not account for the effects on the bubble dynamics themselves of the wall of a tank (Strasberg 1953; Leighton et al. 2002) or pipe (Leighton et al. 1995, 1997, 1998b, 2000; Leighton 2011), or the departure from linearity or steady-state conditions (Clarke & Leighton 2000), all of which could be included given sufficient computational resources (Leighton et al. 2004). Figure 4a plots Yψ, the ratio of the imaginary part to the real part of the complex wavenumber as defined in equation (2.4), for infinite volumes of homogeneous bubbly liquid. The dashed line is the value of Yψ that would have been obtained if the BSD that occurred in the pipe (figure 3, black solid line) were to have existed in just such an infinite body of liquid. The open symbols and asterisks are the result of using equation (2.9) to convert the measured Yψ of the various modes in the pipe into approximate values appropriate for an infinite volume of bubbly liquid. The shape of the ET1 result (triangles) is similar to the dashed curve except an abrupt peak around k1a≈3.4, where the break point (the frequency where the phase speeds of the ET1 and ET2 modes (the solid lines in figure 1a) approach closest with each other, Baik et al. 2010) is observed. In the frequency range of k1a<2, the approximated Yψ is greater than the dashed curve by the factor of about 10 (experienced users will often discard end-of-range data because such effects are expected). Except for these two frequency ranges, the approximated Yψ converted from the ET1 mode is similar to the attenuation in PWFF bubbly liquid. Since equation (2.9) assumes small deviations between axial wavenumbers, q′0mf and q′0mb, and between free field wavenumbers, k′1f and k′1b, wherever these conditions are disobeyed, the approximation inherent in equation (2.9) does not work well. Therefore, conversion from the ET1 mode gives a large discrepancy near k1a≈3.4 where the difference between q′0mf and q′0mb becomes greatest. This can be seen from figure 1, recalling that the normalized phase speed C0m/C1 is calculated from k1/Re[q′0m], where Re represents real part and q′0m is either q′0mf or q′0mb. In figure 1, k1a≈3.4 represents the location where the separation between the ET1 solid and dashed sound speed curves in (a) becomes greatest, and where attenuation in (b) peaks. A similar large discrepancy between these two curves is also observed in figure 1a for k1a<2 where the difference between k′1f and k′1b is large for the ET1 mode. Use of the Yψ data from the other modes (ET2 (squares) and ET3 (asterisks)) produces poorer agreement. The ratio Yψ converted from the ET2 mode (squares) is far from the free field Yψ for an infinite bubbly liquid (dashed curve) and, moreover, it gives negative values for Yψ around k1a∼4.4 (which are not displayed in this picture). The same trend is investigated in the ET3 result around k1a∼8. This level of agreement illustrates the difficulties inherent in comparing the values of Yψ measured in real scenarios with those predicted by PWFF theory.
The converted Yψ plot is then used in the acoustic inversion using linear matrix representation in equation (2.3). The input ratio, Yψ, represents the bubble-induced change of complex wavenumber between in bulk bubbly liquid and in bulk bubble-free liquid (both being in PWFF conditions). For laboratory experiments using the frequency span and source-to-receiver separations used here, the attenuation in an infinite body of bubble-free liquid is usually negligible (at least 1000 times smaller than that in bubbly water), which means Yψ is very small. Since use of equation (2.9) shows greatest success in converting the Yψ of the ET1 mode, it is not surprising that when this converted Yψ (triangles) is used in the inversion (that assumes an infinite body of homogeneous bubbly liquid) in figure 4b, it is most successful at predicting the actual BSD that was measured optically. Indeed, the level of agreement is similar to that shown here for the standard technique, when data from actual PWFF conditions (the dashed curves in figure 1c,d) is inverted using the PWFF inversion to generate the BSD shown by the open circles in figure 3. The other modes (squares for ET2, asterisks for ET3) are, unsurprisingly, less successful in their estimates in figure 4b of the BSD, because their estimates of Yψ were less good in figure 4a.
It is sensible in such inversions to eliminate unphysical peaks (such as that seen in the ET1 mode (triangles) around k1a≈3.4), typical of inversion artefacts, and fill the resulting absent region with the interpolated Yψ using the remaining smooth raw data (greater regularization would eliminate such artefacts but smooths data unnecessarily, so good practice is to add just enough regularization to make the last negative bubble count turn positive, then eliminate likely artefact peaks).
Using identical symbols schemes, figures 3 and 4b, respectively, show the BSDs obtained from the processed (by equation (2.9)) and the unprocessed modes. Of the three modes tested, the ET1 mode (triangles) produces an estimate from this new approach (triangles in figure 4b) that is closer to the actual BSD than simply using the unprocessed ET1 attenuation (triangles in figure 3). While the BSDs obtained from the processed ET2 and the ET3 modes still give poor results, as expected, they illustrate the cause of this error: the fact that in order to obtain a good estimate of BSD, the lowest frequency for which data is available must be lower than the pulsation resonance of the largest bubble present in the population, as discussed earlier. The cut-off frequencies of the ET2 and the ET3 modes are around k1a=3 and k1a=5.5, respectively, corresponding to resonance frequencies of 200 and 100 μm radius bubbles. Because there is therefore no information for Yψ at lower frequencies than these cut-offs for each mode, the acoustic inversion of the modes only gives finite bubble counts in figure 4b for bubbles smaller than the radii corresponding to these cut-off frequencies. Although it is possible to obtain estimates of the BSD above these radii, those are expected to be highly ill-posed despite regularization. Therefore, in order to get the best estimation of the BSD, it is necessary to have a measurement of the Yψ at frequencies less than the resonance frequency of the largest radius of the bubble in the BSD. This is a challenging requirement experimentally, since for a given broadband source, the cost per octave is roughly constant with today's manufacturers, making the low-frequency region particularly problematic: it costs roughly as much to cover the 1.5–3 kHz band as it does to cover the 20–40 kHz band. To save cost, experimenters will choose a low-frequency limit for their sources, but should not then apply this system unless they know that no bubbles with resonances lower than this low-frequency limit will be present. Use of optical techniques to measure the larger bubble sizes, can mitigate this problem. Optical techniques can fail to resolve the smallest bubbles but become more accurate and cost-effective for larger ones. In this way, they exhibit the opposite trend to acoustic methods, and are therefore complementary to acoustic techniques. The proposition that parametric sonar might be applied to explore the low-frequency regime is complicated by the fact that bubbles affect the parametric sonar operation itself, not just the wave once generated.
(b) Acoustic measurements
Section 4a illustrated how difficult is the problem of inferring BSDs using acoustic inversion in a pipe when the data are simulated (and so covers the full frequency range without gaps). Real measurements produce less idealized coverage, and hence the problems of using it in an inversion are greater. The original design for this system commissioned transducers (that would wet when immersed in mercury, but not react with it) covering the 1 kHz–1 MHz frequency range, but this was no longer affordable after the 2008 global financial crash. Therefore, the current exercise was conducted to see if the affordable bandwidth was worth purchasing. Figure 5 shows the measurements of phase speed and the attenuation of the ET1 (triangles), ET2 (squares) and ET3 (asterisks) modes taken in the bubbly water-filled PMMA pipe. The theoretical dispersion curves under the BSD in figure 3 (black solid curve) are calculated from the theory of Baik et al. (2010) and plotted for comparison as solid curves. Although all three modes produce theoretical curves across the frequency range studied, in practice certain modes are observed preferentially in particular frequency ranges for the reasons discussed by Baik et al. (2010). It should be noted that although the pipe used here is the same as that used in Baik et al. (2010), the frequency ranges over which modes are observed change because figure 5 is for bubbly water whereas Baik et al. (2010) made measurements in bubble-free water. The ET1 mode is only observed above the frequency of k1a∼3.5. Below that frequency, the ET2 mode is observed instead. The ET2 mode is again observed above k1a∼6.5. The ET3 mode is observed below this frequency and its phase speed becomes infinite when the acoustic excitation is close to the cut-off frequency of the ET3 mode. Compared with the agreement seen in the ET2 and ET3 measurements, the phase speed of the ET1 modes in figure 5a was measured to be generally less than the theoretical prediction.
Figure 5b shows the measurement of the attenuation of the modes, using the same symbol notation as for figure 5a. As discussed by Baik et al. (2010), the attenuation obtained from Prony's method is sensitive to the sampling range along the axis of the tube (where the z1 and z2 locations of hydrophone along the tube axis represent the start and endpoints of the sampling). The error bars represent ±1s.d. in the mean of results obtained by taking different sampling ranges for |z2−z1| from 0.3 to 1.3 m in 1 cm increments. Although error bars are also plotted in figure 5a, those are very small and mostly extend inside the markers (Baik et al. 2010). Although at first sight, the experimental data appear to follow the theoretical predictions, validating the theory, the degree of validation must be tempered by the large size of the error bars, as was predicted by Baik et al. (2010). In reality, the large error bars imply that all that can conclusively be stated is that the data do not contradict the theory proposed in this study.
Since Yψ=υψ/uψ is the quantity which is used in the inversion, attenuation associated with υψ must be neither so great that there is poor SNR, nor so small that the value of the Yψ parameter fed into the inversion has a proportionally massive uncertainty. Because of this, and because it does not have a cut-off frequency (and therefore incorporates information about the large bubbles), the ET1 mode is the most useful of those measured here for inferring the BSD. In fact, a peculiarity of the actual measurement conducted here relates to the issue mentioned above, whereby unlike the simulated data tests of the previous section, the real dataset is incomplete. The peculiar consequence here is that, in the frequency range over which signal for ET1 is detectable (3.9<k1a<6.2), its attenuation is almost identical to that which the optically measured BSD would give in an infinite body of bubbly water (figure 1). Therefore, one would not expect this ET1 data to be affected much by conversion through equation (2.9), and then to give a reasonable estimate of BSD over a limited range of bubble sizes. Figure 6a shows the converted Yψ from the average of measured Yψ (ET1 (triangles), ET2 (squares) and ET3 (asterisks)) by the process in equation (2.9). The theoretical Yψ in the infinite body of bubbly water with the measured BSD inside the pipe (black solid curve in figure 6b) is represented by solid curve. When these data in figure 6a are used in the acoustic inversion, the result is displayed using the same symbols in figure 6b. Since the measurement of the ET1 mode is only obtained in the frequency range of 3.9<k1a<6.2 (equivalent to 20.6 to 32.8 kHz), the corresponding BSD estimation spans only from approximately 100 to 150 μm radius, which are bubble sizes giving those resonance frequencies. Therefore, the BSD obtained from our measurement is limited to this bubble size range. For such small size ranges, the end-of-range errors that are known to occur with such inversions (and cause divergence of the result at the end of the range) are given undue prominence: as expected therefore, closest agreement between the optically measure BSD and that obtained from the acoustic inversion occurs in the middle of the acoustic range. In order to obtain the BSD of the other radii of the bubbles, it is necessary to measure the Yψ over a wider bandwidth. The estimated BSD obtained from the measurement of the bubbly cylinder is greater than the actual BSD by the factor of from 2 to 4. In this frequency range, the same trend is observed in figure 4b when the theoretical ET1 Yψ in bubbly cylinder was directly used in the acoustic inversion. However, as mentioned above, the theoretical ET1 Yψ in this frequency range is not far from the free field Yψ in bulk bubbly water. When the theoretical free field Yψ was used in the acoustic inversion instead, the estimated BSD is close to the actual BSD as shown by the open circles in figure 3. Therefore, the BSD in this corresponding radius range (from approx. 100 to 150 μm) is very sensitive to the Yψ in the corresponding frequency range (from 20.6 to 32.8 kHz). Consequently, a more precise measurement is required to estimate the correct BSD by the acoustic inversion. However, more precise measurements in this bubble radius range (or corresponding frequency range) do not significantly change the gas volume void fraction since the greatest contribution to the void fraction comes from the larger bubbles, the resonance frequencies of which are far below this frequency range.
The BSDs estimated from the higher modes ET2 (squares) and ET3 (asterisks) are also plotted in figure 6b. One might expect a poor estimate of the BSD from ET3, because the converted values of Yψ for this mode (the asterisks in figure 6a) are far from the theoretical free field Yψ in bulk bubbly water. Although part of the measurement of ET3 is obtained beyond the frequency range of the measured ET1 mode (k1a∼6.4), and so might have been potentially useful in extending the range of bubble sizes estimated, the accuracy of those estimations (the asterisks in figure 6b) do not agree well with the actual BSD. (Indeed, although the BSD estimated from the ET3 mode spans to a higher frequency range than that covered by the ET1 mode, the BSD data points corresponding to these smaller bubbles sizes are not plotted in this figure because (at ) the values are out of range of the plotted bubble concentrations). The same trend is observed for ET2 mode when the data at k1a∼6.6 were used for the acoustic inversion (the corresponding frequency resonates with bubbles of approx. 90 μm radius), which gives BSD estimations (squares) that do not agree with the solid BSD curve. However, the converted values of Yψ from the measurement of the ET2 mode (squares in figure 6a) are very close to the theoretical free field Yψ for a bulk bubbly liquid in the range of 2.9<k1a<3.9 which corresponds to frequencies that are resonant with bubbles having radii from approximately 150 to 210 μm. Consequently, the BSD estimated from the ET2 mode in the bubble radius range of approximately 150–210 μm (the two squares in figure 6b) are close to the black solid BSD curve. As expected, there is much poorer agreement for the estimations of BSD from bubbles that are smaller than approximately 150 μm radius for the ET2 mode. Because of the missing data in the range of 3.9<k1a<6.5, these estimates were obtained from the interpolation using existing ET2 data, which is, of course, far from the theoretical Yψ, and therefore gives poor agreement with the actual BSD over the corresponding bubble size range (radii from approx. 90 to 150 μm).
Simulated and real data were used to explore how measured phase speeds and attenuations in bubbly liquid in a pipe might be inverted to estimate the BSD (which was independently measured using an optical technique). Use of the standard PWFF theory generated large errors when applied to data in pipes. A new inversion used a Taylor series to relate the complex wavenumber measured in a pipe to what would, for the same BSD, have been measured in an infinite homogenous volume of bubbly liquid. When simulated data for the frequency range originally planned for the sensor were used, the new inversion performed well. However, it performed poorly when using actual data from the more narrow frequency range that could be afforded, the result being dominated by the known end-of-range errors that occur during such inversions. The conclusion is that it is not worth ORNL commissioning a sensor with this reduced frequency range, because although it is affordable, the loss of accuracy is unacceptable While such errors can be identified by using complementary optical techniques, this is not possible in optically opaque media and pipelines (although the authors developed an optical fibre which detected bubbles in mercury, a system which provided statistically significant sampling was unaffordable). It is therefore important to question the accuracy of estimated BSDs particularly when they inform high-value decisions by industry. The planned addition of a controlled population of helium bubbles to ORNL's SNS TTF in 2013 is designed to mitigate damage to the target, and so prolong the period between which the target is replaced. The cost of a failed target could amount to £12 M, including $800 k for the new target, $200 k for disposal of the old (now radioactive) target, and the cost of interruption to the use of the facility (B. Riemer 2012, personal communication). The method in this paper allowed industry to assess beforehand the cost/benefit of employing a given inversion technique, including the cost of insufficient investment in a wide enough frequency range.
This work is supported by the ORNL, TN, USA. ORNL is managed by UT-Battelle, LLC, under contract DE-AC05-00OR22725 for the US Department of Energy; and by the UK Science and Technology Facilities Council (STFC). The authors thank B. Riemer and M. Wendel (ORNL) and C. Densham, O. Caretta and T. Davenne (STFC) and Ron Roy of Boston University for advice, discussions and access to facilities.
- Received February 2, 2012.
- Accepted March 29, 2012.
- This journal is © 2012 The Royal Society