Determination of Particle Size Distributions from Acoustic Wave Propagation Measurements

The wave equations for the interior and exterior of the particles are ensemble averaged and combined with an analysis by Allegra and Hawley (cid:64) J. Acoust. Soc. Am. 51 , 1545 (cid:126) 1972 (cid:33)(cid:35) for the interaction of a single particle with the incident wave to determine the phase speed and attenuation of sound waves propagating through dilute slurries. The theory is shown to compare very well with the measured attenuation. The inverse problem, i.e., the problem of determining the particle size distribution given the attenuation as a function of frequency, is examined using regularization techniques that have been successful for bubbly liquids. It is shown that, unlike the bubbly liquids, the success of solving the inverse problem is limited since it depends strongly on the nature of particles and the frequency range used in inverse calculations. © 1999 American Institute of Physics. (cid:64) S1070-6631 99 (cid:33) 01405-1


I. INTRODUCTION
Determining the particle size distribution of a solidliquid mixture is of great practical interest. It has been suggested in the literature that this distribution may be determined by measuring the attenuation of a sound wave propagating through the mixture as a function of the frequency of the wave. The main premise is that the attenuation caused by a particle as a function of frequency depends on its size and therefore the attenuation measurements can be inverted to determine the particle size distribution-at least when the total volume fraction of the solids is small enough so that the particle interactions and detailed microstructure of the slurry play an insignificant role in determining the acoustic response of the slurry. Indeed, this general principle has been exploited successfully to determine the size distribution of bubbles in bubbly liquids. 1-3 Commercial ''particle sizers'' based on acoustic response are in the process of being developed/marketed for characterizing solid-liquid mixtures. 4 The main objective of this paper is to investigate under what circumstances such a problem can be solved for solid-liquid systems. It will be shown that the success of the acoustic method for determining detailed particle size distributions is limited, depending on the nature of the particles and the frequency range over which input data ͑attenuation͒ are available.
The problem of determining the acoustic response of a slurry given its particle size distribution is referred to as the forward problem. When the total volume fraction of the particles is small, the problem is relatively simple since then one only needs to understand the interaction between a single particle and the incident sound wave. This has been examined by a number of investigators in the past with notable contributions from Allegra and Hawley 5 and Epstein and Carhart 6 who considered suspensions of particles as well as drops. The former investigators also reported experimental results verifying the theory for relatively small particles for which the acoustic wavelength is large compared with the particle radius. The theory developed by these investigators is quite general and accounts for attenuation by thermal, viscous, and scattering effects as described in more detail in Secs. II and III. The case of monodisperse nondilute suspensions has been examined by Varadan et al. 7 who used an effective medium approximation to account for particle interactions, but their analysis was concerned only with the attenuation due to scattering. In Sec. II we present the theory for the forward problem with the main aim of reviewing the important physical effects causing the attenuation. The derivation for the attenuation proceeds along different lines than that followed by Epstein and Carhart or Allegra and Hawley in the way the one particle analysis is used to predict the attenuation of the suspension. These investigators calculated the energy dissipation per one wavelength to estimate the attenuation while we use the method of ensemble averages to determine both the phase speed and attenuation of waves. The method of ensemble averages will be more convenient for developing a suitable expression for attenuation in nondilute suspensions, if desired, using either an appropriate effective-medium approximation or direct numerical simulations.
In Sec. III we present new experimental data for nearly monodisperse polystyrene particles whose radii are comparable to the wavelength and validate the theory described in Sec. II over a nondimensional frequency range much broader than examined by previous investigators. We also summarize in that section the different physical mechanisms that cause attenuation in suspensions. The attenuation as a function of frequency is shown to undergo several peaks owing to the resonances in shape oscillations in agreement with the theory prediction. It also gives some indication of the range of frequency and attenuation measurable with our acoustic device.
In Sec. IV we consider the inverse problem, i.e., the problem of determining the particle size distribution given the total attenuation as a function of frequency and the physical properties of the particles and the suspending liquid. At small particle volume fractions, this amounts to solving a linear integral equation for the unknown size distribution. This is an ill-posed problem: small changes/errors in the attenuation data can cause large changes in the size distribution. Thus, for example, several very different particle distributions could give rise to essentially the same attenuationfrequency curve. This, of course, is a rather well-known difficulty in most inverse problems which involve solving a Fredholm integral equation of the first kind with a smooth kernel. Techniques have been developed to ''regularize'' the problem. We use the well-known Tikhonov regularization techniques, 8 which replaces the ill-posed original problem with another well-posed problem involving an integrodifferential equation whose solution minimizes the fluctuations in the particle size distribution. Minimizing of the fluctuations is rationalized on the grounds that in most practical situations the particle size distribution is smooth. This regularization technique has been shown to work very well for the inverse problem in bubbly liquids. 2 We apply the above technique to suspensions of polystyrene and glass particles. We find that the technique works well for the polystyrene particles but not for all glass particles. We also find that for polystyrene particles the technique works only when the attenuation is given over an appropriate frequency range-a frequency range that is too narrow or too broad may give erroneous estimates of the distribution. An alternative inverse technique based on linear programing also failed to produce the correct particle size distribution for the cases for which the Tikhonov scheme failed. This suggests that the prospects for determining the detailed particle size distribution from acoustic measurements are somewhat limited. ͑In situations where more might be known about the nature of particle size distribution, e.g., unimodal with a Gaussian or log-normal distribution, one might be able to determine the size distribution through appropriate curve fitting as has been done, for example, by McClements and Coupland, 9 but this is not addressed here.͒ The reasons why the size distributions for some particle suspensions are not recovered by the inverse techniques while the same techniques were found to be quite successful for bubble suspensions can be given in terms of differing resonance nature of these suspensions. In the case of bubbles in most typical applications, the resonance occurs at frequencies for which the wavelength is relatively large compared with the bubble radius. This resonance is due to volume oscillations; the shape-dependent resonances are unimportant and, as a consequence, there is effectively one resonance frequency for each bubble size. Thus, the peaks in the attenuation-frequency curve give a reasonable indication of the bubble sizes. The situation with the particles is different as their resonance behavior is governed by shape oscillations. For polystyrene particles, several resonance peaks corresponding to different shape oscillations arise even for monodisperse particles, and, as a result, it is difficult to determine whether a given resonance peak arises from a different shape oscillation mode of the same particle or from a particle of different size. For glass particles, on the other hand, there are no significant resonance peaks even for monodisperse particles, and the attenuation behavior for different sizes is not significantly different to allow accurate results for the size distribution.

II. THE FORWARD PROBLEM
The wave equations for both the interior and exterior of particles have been derived by Epstein and Carhart. 6 They were interested in the attenuation of sound waves in fog and therefore their analysis was concerned with drops instead of particles. The stress tensor for a viscous fluid used by them for the interior of the drops was subsequently replaced by Allegra and Hawley 5 by that of an elastic solid to determine the attenuation of sound waves in a solid-liquid suspension. In this section we shall ensemble average a resulting wave equation to obtain the effective wave number of the suspension at arbitrary volume fraction, the real and imaginary parts of which give the wave speed and attenuation. Thus, the attenuation is not calculated by means of an energy dissipation argument, 5,6 but directly from averaging the relevant wave equation. The result contains certain coefficients that remain to be evaluated for a given microstructure. In the present study, since we are primarily concerned with determining the size distribution, we shall evaluate the coefficients in the limit of small volume fractions. In a separate study, where we shall present experimental results for nondilute suspensions, we shall extend the theory to treat nondilute suspensions.

A. Theory
Epstein and Carhart 6 first linearized the conservation equations for mass, momentum, and energy. The pressure and internal energy were then eliminated by introducing the linearized equations of state to yield equations in terms of density, velocity, and temperature. Next, the time dependence of all quantities were expressed by the factor exp(Ϫt)-which is henceforth suppressed-and the velocity was expressed as vϭϪٌ⌽ϩٌϫA, with ٌ•Aϭ0. With this form of v it is possible to eliminate the temperature and density from the governing equations to yield a fourth-order partial differential equation for ⌽ and a second-order equation in A. The former, in turn, can be split into two second-order wave equations upon a substitution ⌽ϭ c ϩ T to finally yield three wave equations: The wave numbers in the above equations are given by Here, c is the phase speed in pure liquid, is the density, and are, respectively, the compressional and dynamic viscosity, ␥ϭC p /C v is the ratio of specific heats at constant pressure and volume, is the thermal conductivity, and ϭ/C p is the thermal diffusivity. Inside the particles similar equations hold with the dynamic viscosity replaced by /(Ϫ) and the wave speed by "( ϩ2 /3)/ … 1/2 , where and are the Lamé constants, while the compressional viscosity is left out. Henceforth a tilde refers to the interior of particles.
At small values of e and f ͑such as in water͒, the above expressions for k c and k T simplify to and its counterpart inside the particles describe the sound wave propagation through the suspension. Note that the wave number has an imaginary part; sound waves in pure fluid are attenuated by viscous and thermal energy dissipation; 10 the term inside the square brackets in ͑8͒ is commonly referred to as the ''diffusivity of sound.'' The total attenuation coefficient in both liquid and in the solid particle will henceforth be treated as additional physical properties. The other two wave equations describe waves that arise from thermal conduction and finite viscosity: we note that the modulus of k T in Eq. ͑8͒ is inversely proportional to the thermal penetration depth ͱ/ and that of k s to the viscous penetration depth ͱ/. The thermal ( T ) and shear ͑A͒ waves have generally very high attenuation and are unimportant in acoustic applications.
We now proceed to ensemble average the wave equation ͑1͒ to find an expression for the effective wave number of a wave propagating through a solid-liquid suspension. Introducing an indicator function g, defined to be unity inside the particles and 0 outside, the ensemble-averaged value of c is To obtain a wave equation for ͗ c ͘ we first take the gradient of the above equation to yield ٌ͗ c ͘ϭ͗gٌ c ϩ͑1Ϫg ٌ͒ c ͘ϩٌ͗͑g͒͑ c Ϫ c ͒͘. ͑9͒ As argued by Sangani,11 upon assuming that the particles' spatial distribution is homogeneous on a macroscale, the last term in ͑9͒, being a vector, can only depend on quantities such as ٌ͗ c ͘ and ٌٌ 2 ͗ c ͘. Anticipating that ͗ c ͘ will satisfy a wave equation we express the last term on the righthand side of the above equation in terms of ٌ͗ c ͘, i.e., we write ٌ͗͑g͒͑ c Ϫ c ͒͘ϭ 1 ٌ͗ c ͘, where 1 depends on the parameters such as the volume fraction, k c , and k c . The divergence of ͑9͒ is given by we find that ͗ c ͘ satisfies a wave equation with the effective wave number given by The real part of the effective wave number is the frequency divided by the phase speed in the mixture and the imaginary part the attenuation. Up to this point the analysis is rigorous and without any assumption. Applying the boundary conditions of continuity of temperature, flux, velocity, and traction at the surface of the particles, and solving the resulting boundary value problem numerically, it is possible, in principle, to determine the phase speed and attenuation at arbitrary volume fraction using the above formulation. Special simplifications can be made when the wavelength is large compared with the particles and when the viscous and thermal depths are small compared with the particle radius for which numerical computations using the multipole expansions developed in recent years ͑see, e.g., Ref. 12͒ can be readily used for determining the attenuation at arbitrary volume fractions. Alternatively, one may use a suitable effective-medium approximation to account for the particle interactions in nondilute suspensions using the above framework. We shall pursue this further in a separate study 13 devoted to nondilute suspensions where we shall also present experimental data. Since our interest in the present study is in determining size distributions, it is necessary to consider only the simplest case of dilute suspensions.
In dilute suspensions the particle interactions can be neglected, and the coefficients 1 -3 can be evaluated from the solution for c for a single particle given by Allegra and Hawley. 5 Accordingly, the conditionally averaged ͗ c ͘(x͉x 1 ) given a particle centered at x 1 is given by where rϭ͉xϪx 1 ͉, ϭcos , being the angle between x Ϫx 1 and k c , h n is the spherical Bessel function of the third kind ͑corresponding to an outgoing scattered wave͒, and P n is the Legendre polynomial of degree n. The first term on the right-hand side of the above expression is the unconditionally averaged ͗ c ͘(x) whose amplitude is taken to be unity with no loss of generality. Inside the particle centered at x 1 we have where j n is the spherical Bessel function of the first kind. Similar expressions are written for the conditionally averaged T and A. This results in expressions with a set of six unknowns for each mode n. Application of the aforementioned boundary conditions of continuity of velocity, traction, temperature, and heat flux yield six equations in these six unknowns for each n. There were some typographical errors in the equations given by Epstein and Carhart 6 and Allegra and Hawley; 5 the correct equations are given in the Appendix. Although it is possible to solve for the unknowns analytically in certain limiting cases, it is best to solve them numerically since we are interested in covering a wide frequency range for inverse calculations.
We now return to the calculations of the coefficients 1 -3 . Upon using the identity ٌgϭϪn␦͑xϪx i ͒, with x i being a point on solid-liquid interface and n the unit normal vector at the point, 1 is given by Here, P(x 1 ) is the probability density for finding a particle in the vicinity of x 1 . Similarly, we have for 2 and 3

͑17͒
The above integrals must be evaluated while keeping in mind that the integration variable is x 1 . Thus, for example, in ͑15͒ and ͑16͒ we must consider all particles whose surfaces pass through the point x. To carry out these integrals we use the identity and the orthogonality of the Legendre polynomials over spherical surfaces. The resulting expressions are where in the expression for 1 the j nϪ1 -term in the nϭ0 contribution is to be left out. Here, is the volume fraction of the solids, zϵk c a and zϵk c a are the nondimensional wavenumbers, and primes denote derivatives. Expressions ͑19͒-͑21͒, together with the expression for the effective wave number ͑12͒, complete the description of a solid-liquid mixture at low volume fractions.
In the above we have derived expressions for the attenuation and wave speed by calculating the effective wave number directly. An alternative derivation of the attenuation coefficient is to calculate the energy dissipation per wavelength in the mixture and divide the result by the energy per wavelength. The result for the attenuation per unit length is then 5,6 ␣ϭϪ 3 It can be shown that the two methods give essentially the same result for the attenuation in the limit →0 with z Ϫ2 Re A n in the above replaced by Re(A n /z)/Re(z) in the ensemble-averaging method presented here. The above analysis may be extended to account for the effect of finite volume fraction through a suitable effectivemedium approximation. Sangani 11 showed that the first correction of O( 3/2 ) to the dilute O() approximation for bubbly liquids can be simply derived through an effectivemedium approximation. This correction is most significant near the resonance frequency of bubbles, and to correctly capture the behavior near resonance it is important to replace the pure liquid wave number ͑k c in the above analysis͒ by the effective wave number. Thus, in the present context, zϵk c a in ͑19͒-͑21͒ for 1 -3 , is replaced by z eff ϵk eff a, while the wave number in pure liquid in the expression for k eff , ͑12͒, has to be retained. The latter expression is then iterated to obtain a converged solution for k eff . The effective-medium approximations have been found to be quite useful in the related study of light scattering by suspensions ͑see, e.g., Ref. 14͒. For very high volume fractions the other physical properties of the so-called effective medium must also be modified. In a separate study, 13 where we shall report experimental data for dense slurries, we shall examine several different versions of effective-medium approximations in more detail.
Finally, the above analysis can be extended in a straightforward manner to account for the particle size distribution when the total volume fraction of the particles is small. Let us write the attenuation by particles of radius between a and aϩda as an attenuation density ␣ ( f ,a) ͓where f is the frequency of the wave, f ϭ/(2)͔ times the volume fraction of those particles (a)da; we shall refer to (a) as the volume fraction distribution. At low volume fractions these contributions can be ''summed'' over all particle sizes to give for the total attenuation ␣ tot (f ): It is customary to express the particle size distribution in terms of its number density distribution P(a). The volume fraction distribution is related to P(a) by (a) ϭ(4a 3 /3) P(a).
The effective-medium approach described earlier can also be readily extended to account for the particle size distribution. The coefficients 1 -3 are first determined as functions of a for an assumed value of the effective wave number and these are integrated after multiplying by (a)da to yield estimates for the average values of 1 -3 for the suspension. These are substituted in ͑12͒ to determine k eff . If this estimate of k eff is different from the the assumed value, then 1 -3 are estimated for the new value of k eff , and the process is repeated until the assumed and estimated values of the effective wave numbers agree with each other.

III. DISCUSSION AND COMPARISON WITH EXPERIMENTAL DATA
Figures 1 and 2 show the predictions for the attenuation and wave speed as a function of frequency f for 79 m radius polystyrene particles at a volume fraction of 0.05. The frequency f in Hz is related to by ϭ2 f . The physical properties used in the computations are given in Table I. 15 We note that the wave speed only changes if the frequency becomes very large and that these changes coincide with strong changes in the attenuation as well. Hence we expect that the measurement of the phase speed will not provide significantly new information over that obtained from the attenuation measurements alone as far as the problem of determining the size distribution is concerned. On the other hand, since the phase speed at low frequencies is nearly independent of the frequency or k c a, it might be possible to use the low frequency speed data to determine the total volume fraction of the particles regardless of its size distribution. We shall focus in the present study on the results for attenuation as they are the most sensitive to the particle size distribution.
The attenuation of sound waves in a suspension is different from that in pure liquid because of four effects. First, the attenuation of sound in pure solid is different from that in pure liquid, and hence simply the presence of the particles changes the attenuation from that of pure liquid. Second, changes in temperature are different in a solid than in a liquid, and this causes a heat flux through the surface of the particles. This heat flux is out of phase with the sound wave passage and this leads to attenuation referred to as the thermal attenuation. Third is the viscous energy dissipation caused due to the motion of the boundary of the suspended particles. Finally, the fourth effect is the attenuation due to scattering.
Allegra and Hawley 5 showed that when the particle size is much smaller than the wavelength and much greater than the thermal and viscous penetration depths (/) 1/2 and (/) 1/2 , the resulting viscous and thermal attenuations increase as f 1/2 . On the other hand, when the penetration depths are much greater than the particles, both attenuation contributions increase as f 2 . This transition occurs at very low frequencies-about 2 Hz for 100 radius particles in water-and will not be considered here. Attenuation due to scattering becomes important when the nondimensional wave number zϭk c a becomes comparable to unity. For small but finite z the scattering losses increase as f 4 . Thus, one expects that the change in the attenuation behavior from f 1/2 at low frequencies to f 4 at high frequencies will provide an important indication of the particle size. These asymptotic ranges are indicated in Fig. 1. We see that the transition to the f 4 behavior does not fully occur for the particles considered here. As the frequency is increased particles undergo several resonances as described in more detail below, and this is responsible for the several peaks seen in Fig. 1. Figure 3 shows the contributions to the total attenuation from each P n mode. The nϭ0 mode corresponds to radial ͑volume͒ oscillations of the particles, nϭ1 to the translational oscillations, nϭ2 to the ellipsoidal P 2 -shape deforma-tion oscillations, and so on. The density of polystyrene particles is essentially the same as that of water, hence the particles' translational oscillations are very small. As a consequence, the viscous attenuation is small and the small frequency behavior is governed by the thermal attenuation of the nϭ0 mode. At higher frequencies the nϭ0 mode begins to increase first as f 4 due to scattering losses but the contribution from the nϭ2 mode soon becomes important as it undergoes a resonance at about 3 MHz frequency. The nϭ3 and nϭ1 modes undergo resonances next, and so on. We see that the nϭ0 mode undergoes a broad maximum around 9 MHz. Although not shown here, it too undergoes a resonance with a sharp downward peak at about 21 MHz. Thus, we see that the attenuation varies with frequency in a rather complicated manner at high frequencies owing to various resonances. We should note here that the behavior of this kind for polystyrene particles has also been reported by other investigators in the past. For example, Anson and Chivers 16 and Ma, Varadan, and Varadan, 14 who restricted their analysis to scattering losses only, found essentially the same behavior, and in earlier investigations 17,18 mainly focusing on the determination of waves reflected by immersed objects, high-amplitude reflected waves were found at certain resonance frequencies. Figure 4 shows attenuation as a function of nondimensional wave number k c a for particles of radii 50 and 79 microns. We see that the curves for these two radii are essentially the same, indicating that, at least for polystyrene particles, the thermal or viscous effects have negligible influence on the resonance frequency. The first resonance corresponding to nϭ2 appears at k c aӍ1. 4.
Allegra and Hawley 5 tested ͑22͒ extensively against their experiments and found very good agreement. However, their particles were always smaller than 1 m radius, so that the wavelength was always much greater than the particle size. No resonance behavior was observed in their experiments. Although the above-mentioned paper by Ma, Varadan, and Varadan 14 presents experimental data on light scattering in the small-wavelength regime, no data on attenuation of sound waves by particles were presented. To test how well the theory works for larger particle sizes, we carried out an experiment that will be described in detail ͑along with more experiments on concentrated slurries͒ elsewhere. 13 In this experiment the attenuation of sound waves was measured in a frequency range of 1-10 MHz in a solid-liquid mixture of polystyrene particles with 79Ϯ3 mean radius and 1.8 standard deviation at 0.05 volume fraction. Monochromatic tonebursts, at incremental frequencies, were transmitted by a transducer on one side of a small vessel in which the mixture was being stirred; a second transducer received the signal and sent it to a LeCroy 9310A digital oscilloscope. The amplitude of the signal for pure water was measured, as was that for the solid-liquid mixture. The excess attenuation was determined by where d is the distance between the transducers and V mix and V H 2 O are the voltage amplitudes of the received signals in the mixture and pure water, respectively. The distance between the transducers was 2 in. at low frequencies and 1 in. at higher frequencies; this was necessary because the attenuation at higher frequencies was too large to produce significant signal-to-noise ratio in the larger vessel. Figure 5 shows the comparison between theory and experiment. At the two gaps in the frequency domain ͑where the theory predicts very high peaks͒ the attenuation became again so large that the signal-to-noise ratio vanished even in the smallest vessel. Good agreement is found between experiments and the theory except near resonance frequencies where small differences are seen. There are two possible reasons for these small differences. The first is concerned with the finite volume-fraction effect. To investigate this we have also plotted in Fig. 5 a result from an effective-medium approach described in the previous section. The resulting attenuation changes, but in the wrong direction. The second reason is that the particles were not exactly monodispersed. Using the method described in the previous section, a log-normal particle size distribution was introduced with a mean radius of 77 and 2.5 m standard deviation, which lies within the manufacturers' specifications. The result for the attenuation, the dashed curve in Fig. 5, shows close agreement with the data. Thus, we conclude that the agreement between the theory and experiment is excellent, and that the small observed differences are due to small polydispersity of the suspension.
The attenuation behavior displayed by polystyrene particles is not generic, as can be seen from Fig. 6 which shows the attenuation behavior for glass particles. Since the density of the glass particles is significantly different from that of water, the glass particles execute significant translational oscillations. As a consequence, the low-frequency behavior is completely governed by the viscous effects and the nϭ1 mode. Note that the small frequency attenuation is about two orders of magnitude greater for glass particles than for the polystyrene particles. Also we see a considerably different behavior at higher frequencies. The attenuation does not seem to peak at several frequencies. Rather, for each n we see broad ''hills'' separated by narrow ''valleys.'' The total attenuation does not appear to go through several resonances. The difference in the behavior for the glass and polystyrene particles at these high frequencies seems to arise mainly from the different elastic properties of the two materials.

IV. THE INVERSE PROBLEM
We now consider the inverse problem: given the total attenuation ␣ tot as a function of f we wish to determine (a) using ͑23͒. The straightforward method of solving the inte- -., theory for monodispersed particles with effective medium correction for finite volume fraction effects; ---, theoretical result with a particle size distribution with a mean particle radius of 77 m and standard deviation of 2.5 m ͑this is within the range specified by the manufacturer͒. gral equations, i.e., discretizing the integral domain into a number of elements and converting the integral equation into a system of linear equations in unknowns (a k ) at a selected number of points a k in the domain, cannot be used since the resulting equations will be ill conditioned. Figure 7 illustrates the ill-posed nature of the problem. Figure 7͑a͒ shows two very different particle distributions whose attenuation versus frequency curves are seen in Fig. 7͑b͒ to be essentially the same. These curves were obtained by starting with a smooth, log-normal particle size distribution ͓the dashed curve in Fig. 7͑a͔͒ and generating the attenuation versus frequency data using the forward theory ͓the circles in Fig.  7͑b͔͒. A 1% random noise was then added to the data and ͑25͒ with ⑀ϭ0, which is equivalent to the integral equation ͑23͒, was subsequently solved to yield the particle size distribution indicated by the solid line in Fig. 7͑a͒. The pluses in Fig. 7͑b͒ correspond to the attenuation determined from the forward theory using the new particle distribution. Note that the attenuation is evaluated with a smaller frequency increment than the one used for the original distribution. We see that the attenuation from the two distributions agree with each other to within 1% for the frequencies marked by circles. The highly oscillatory particle distribution does show an oscillatory behavior in between the frequency increments, particularly at 10 MHz, but these oscillations occur only for a very narrow frequency range and would have been missed altogether had the attenuation been determined only at the input frequencies.

A. Method
Since the true particle distribution is expected to be smooth, we must only allow solutions that are reasonably smooth. There are several ways of accomplishing this. In the present study, we shall use primarily a regularization technique due to Tikhonov 8 which was successfully used for bubbly liquids by Duraiswami. 2 An alternative method is presented at the end of this section. Accordingly, we multiply ͑23͒ with ␣ ( f ,a)d f and integrate over the frequency range to obtain a simpler integral equation in which the right-hand side is only a function of a: where (a min ,a max ) and ( f min ,f max ) are the radius and frequency ranges. The above integral equation is now regular- ized as explained below by adding a small term ⑀( Ϫl 2 Љ) ͑where primes denote derivatives͒ to its left-hand side. Thus, we obtain

K͑a,aЈ͒͑aЈ͒daЈϭb͑a͒, ͑25͒
where l is a suitably chosen lengthscale and K(a,aЈ) is a kernel defined by is an integro-differential equation and needs two boundary conditions. Usual practice is to take the derivative of (a) to be zero at the two end points: Ј͑a min ͒ϭЈ͑ a max ͒ϭ0. ͑27͒ Note that a min and a max are not known a priori in general. One expects to be zero also at the two end points. Thus, the range (a min Ϫa max ) must be determined by trial and error so that both and its derivatives are approximately zero at the extreme values of a. Now it can be shown that the solution of ͑25͒ subject to the boundary conditions given by ͑27͒ minimizes where E is the measure of error between the actual attenuation and the computed attenuation: Since both E and the second term in ͑28͒, i.e., the integral, are non-negative, minimization of ͑28͒ ensures that the solution of ͑25͒ will be free from large oscillations in . In other words, highly oscillatory distributions such as the one shown in Fig. 7͑a͒ are rendered inadmissible when ͑25͒ is solved with finite, positive ⑀ in place of the original integral equation ͑24͒. Thus, we have regularized the problem of determining .
If we choose a large ⑀, then we decrease the oscillations in but increase the error in (a) since then the equation solved is significantly different from the original integral equation. Small ⑀, on the other hand, yields unrealistic (a) having large oscillations when the data ␣ tot (f ) are not exact. An optimum choice of ⑀ then depends on the magnitude of uncertainty/error in the attenuation-frequency data. In the calculations we shall present here the exact ␣ tot (f ) is first determined using the forward theory for a given (a) and a small random noise of about 1% magnitude is added to it before the inverse calculations are made ͑the effect of noise magnitude is discussed below͒. Thus, we have an estimate of the error in the data, but in general this estimate may not be known reasonably accurately. To determine the optimum ⑀, we solve ͑25͒ for several different ⑀'s and plot E versus ⑀ to find a minimum in E. This, however, may lead to distributions in which (a) may have unphysical negative values for some a. The constraint (a)у0 for all a is satisfied a posteriori by setting (a)ϭ0 for all a's for which the solu-tion of ͑25͒ gave negative values of . The computed value of E for a given ⑀ is then based on (a)у0.
The integro-differential equation ͑25͒ was solved as follows. After discretizing the domain (a min Ϫa max ) into NϪ1 equal segments and the frequency domain into M Ϫ1 logarithmically equal segments we first evaluate the kernel K(a i ,a j ) for i, jϭ1,2,...,N ͓cf. ͑26͔͒ using a trapezoidal rule for the integration over the frequency range. As pointed out by Duraiswami, 2 it is essential to calculate the integral over particle radius very accurately. We assume that (a) varied in a piecewise linear manner in each segment and use a 12point Gauss-Legendre quadrature to evaluate the integral in ͑25͒. A second-order central difference formula was used to evaluate Љ(a) at all points except the end points a min and a max . The boundary conditions Ј(a min )ϭ0 and Ј(a max ) ϭ0 were approximated using, respectively, second-order forward and backward difference formulas. Application of ͑25͒ at all the discretization points together with the boundary conditions can be expressed with the above scheme as a system of linear equations: where j ϭ(a j ) and b i ϭb(a i ). The above set of equations was normalized by dividing all the equations with the greatest element of the kernel K(a i ,a j ), K m for all i,j, times the segment length ⌬aϭ(a max Ϫa min )/(NϪ1). This set of equations was subsequently solved using a standard IMSL subroutine for linear equations. Once j are determined for a selected value of ⑀, we satisfy the constraint j у0 by setting, as mentioned earlier, j ϭ0 for all negative j . The error E as given by ͑29͒ was subsequently evaluated using a trapezoidal rule for integration over the frequency range. The optimum value of ⑀ was determined by stepping logarithmically through several values of ⑀ and plotting E versus ⑀.
A typical result ͑Nϭ30, M ϭ112, f min ϭ0.1 MHz, f max ϭ17 MHz, a min ϭ15 m and a max ϭ35 m͒ for the error E in the resulting attenuation as a function of ⑀ is shown in Fig. 8. Note that ⑀ here is the actual ⑀ divided by K m ⌬a. We see a clearly defined optimum value of ⑀. Computations were also made with larger M to confirm that the resulting volume fraction distribution was not affected by the further refinement in the integration over the frequency range. A remark should also be made of the choice for the length l in ͑25͒. We may regard both ⑀ and l as parameters to be chosen so as to minimize the error E. Taking lϭ(a max Ϫa min )/n we computed E by varying both ⑀ and n with n varied from 1 to N. The three-dimensional plot of E versus n and ⑀ showed that E was much more sensitive to the choice of ⑀ than it was to n. In general, the results with n close to N were slightly better than with those near nϭ1. Based on this observation we chose nϭ30. For larger values of N(NϾ40) we found that choosing nϭN led to more oscillatory behavior for j . This is to be expected since choosing larger n, and, hence, smaller l, permits larger values of Ј(a).

B. Results and discussion
We now present results for the volume fraction distribution obtained using the above technique. As mentioned earlier, we used the forward theory to generate attenuation data for an assumed volume fraction distribution. Small random noise can be added to the data thus generated to mimic possible errors arising in the attenuation measurement. This is satisfactory since we are primarily interested in assessing the procedure for solving the inverse problem. If the procedure gives erroneous results even for this case, it will certainly break down in practice using the experimentally generated data.
The frequency range over which the attenuation measurements are carried out in our laboratory is 0.1-15 MHz. We shall choose here the same range to investigate the success and limitations of the above technique to solve the inverse problem although we shall also consider cases with a larger frequency range to inquire if better estimates of (a) could be achieved if the attenuation data at higher frequencies were to be made available. This is important since the acoustic instruments operating up to 150 MHz are available.
We consider first particle sizes that are of the same order of magnitude as the wavelength somewhere in this frequency range, which is the case for particles of about 10-100 radius ͑for larger particles observed behavior of the attenuation is shifted to lower frequencies͒. A particle size distribution that is often used is a log-normal distribution, which results in volume fraction distributions such as the smooth one shown in Fig. 7͑a͒. We attempt therefore to recover that distribution from the corresponding attenuation data. As in the forward problem, we shall investigate polystyrene particles and glass particles in water, as the first are almost neutrally buoyant and deformable while the latter are very rigid and much heavier than water; the physical properties used in the present calculations are listed in Table I.
We begin with the results for polystyrene particles with a narrow size distribution in the range of 20-30 m. The particle size range for the inverse calculations is first taken to be much greater-5-100 m; the frequency range was 0.1-17 MHz. Figure 9 shows that the volume fraction distribution as evaluated from the inverse technique is in very good agreement with the input distribution. The result for the size distribution can be improved further by making the particle size range smaller ͑a close-up of the improved result is shown in Fig. 11͒.
In Fig. 10 we consider a more complicated, bimodal size distribution in the range of 20-45 m with peaks at about 25 and 38 m. The attenuation as a function of frequency for this distribution is shown in Fig. 10a. The maximum frequency used for inverse calculations is indicated by a square; it is seen that the frequency range includes the first two resonance peaks of the attenuation curve. From Fig. 10͑b͒ we see once again that the inverse procedure recovers this distribution very well.
One of the difficulties in solving an ill-posed problem is that small errors in the input ͑attenuation͒ data can cause large changes in the solution. Of course, errors are always present in the experimentally obtained attenuation data. The calculations presented so far were made with no added noise. To mimic the practical situation, we added random noise of 5% standard deviation to the input data; this is about the same as the order of magnitude of the errors present in the experimental results shown in Fig. 5. The resulting volume fraction distribution, shown in Fig. 11, does confirm that small fluctuations in the input data only cause small deviations in the output. When the calculations were repeated with FIG. 8. Typical dependence of the error in the attenuation for the solved volume fraction distribution as a function of the regularization parameter ⑀. The ͑small͒ parameter ⑀ should be chosen such that this error is minimized. The minimum was always found to be well-defined.
FIG. 9. Solving the inverse problem for polystyrene particles. The solid line is the volume fraction distribution used to generate attenuation data ͓shown in Fig. 12͑a͒, with f max as indicated by a square͔; the dashed line is the solution of the inverse problem when taking the particle radius range to be 1-100 m and using 50 ''bins'' of particle sizes. a noise of 10% standard deviation, the computed particle size distribution was found to be considerably different from the input distribution, although the main features of the size distribution were preserved by the inverse computations.
The results discussed so far suggest that the inverse problem can be solved with reasonable success. We now illustrate some limitations. The inverse method gave erroneous particle size distributions for smaller particles when the same frequency range as the above was used. Of course, in order that the size of the particles be determined there must be at least one transition-from the thermal attenuation dominated regime to the scattering dominated regime which occurs roughly speaking at k c aϭO(1). If the particles are very small, then this transition may not occur over a fixed frequency range. However, as we shall presently see, the results are very sensitive to the frequency range chosen for computations even when this transition is included in the range. Figure 12 shows the effect of varying f max on the computed distribution. As seen in the figure the resonance in the shape oscillations of the ͑polystyrene͒ particles leads to a change in the slope of the curve just before the first resonance. This transition occurs just beyond the point marked by a circle in Fig. 12͑a͒. We see a marked improvement in the results in Fig. 12͑b͒ when f max is chosen corresponding to a point marked plus in Fig. 12͑a͒ over those obtained with a point corresponding to the circle which does not include the second change in slope. The point marked plus corresponds to a frequency greater than the frequency at which the second change in slope occurs for larger particles but smaller than that for smaller particles. This seems to give rise to an inverse solution which is reasonably accurate for larger particles but not for smaller particles. Also shown in Fig. 12͑b͒ are the results when f max is chosen to coincide with the end of first peak, the point marked square in Fig. 12͑a͒. This is seen to yield very accurate results for the size distribution.
One might suppose that covering a broad enough frequency range will alleviate the difficulties seen above. This, unfortunately, is not the case. Figure 13 shows the results for three different f max . The dashed curve corresponds to cutting off the frequency range at the end of first peak as in Fig. 12, the dashed-dotted line to the end of three peaks, and the dotted line to 10 9 Hz, a frequency about 50 times greater than the first resonance frequency. We see that the results of inverse calculations actually deteriorate if a much larger range of frequency is employed, notwithstanding the fact that measurements over such a broad frequency range could itself be a very challenging task. One may rationalize this result as follows. As seen in Fig. 1 a monodisperse suspension will exhibit several resonance frequencies corresponding to various shape oscillation P n (nϭ2,3,...) modes. Thus, a second peak in the attenuation-frequency curve for polystyrene particles may correspond either to say, a P 3 mode of a larger particle, or may correspond to a P 2 mode of a smaller particle. In our calculations we used only up to the first six modes (nр5), but in practice the acoustic response may be further complicated by the higher-order modes for frequencies of order 10 9 Hz considered here.
Since including a wide frequency range with several resonance peaks seems to adversely affect the inverse calculation, one may consider cutting off the attenuation data beyond first peak. This, however, may not work if the distribution is truly bimodal as was the case considered earlier in Fig. 10. If we omit the second resonance peak from the attenuation data by considering a maximum frequency that is less than the point marked square in Fig. 10͑a͒, say, that marked by the circle, we get a poor inversion as shown in Fig. 14. The inverse technique computes accurately the volume fraction distribution of larger particles whose resonance was included in the data but fails to predict that for smaller particles. Figure 15 shows results for a broad, unimodal distribution. The resonance peaks of different particles overlap in this case resulting in the absence of peaks in the attenuationfrequency curve ͓Fig. 15͑a͔͒. Figure 15͑b͒ shows the results of inversion for three different cut-off frequencies. The largest frequency, marked by a square in Fig. 15͑a͒, is larger than the second transition frequency of small as well as large particles, and this seems to produce excellent inverse results.
In most of the inverse calculations shown so far which yielded poor results, we note that the failure is particularly severe for smaller particles. One may rationalize this by observing that the total error E will be dominated by the errors at large frequencies since the attenuation there is very large.
FIG. 12. Influence of the size of the frequency range over which attenuation is specified on the solution of the inverse problem. Polystyrene particles. ͑a͒ Input-attenuation data and four different upper bounds on the frequency. ͑b͒ Results from the inverse problem from these different ranges, using the same marker type. The solid line is the exact result; ᮀ, result when cutting off the frequency range just at the end of the first peak in the attenuation; ϩ, result when cutting of the frequency range after the second change in slope of the attenuation; and ᭺, result when cutting off before the second change in slope. Cutting off the frequency range at the point marked ''छ'' is discussed along with Fig. 13. When k c a min Ͻ1 in the frequency domain that is considered, the small particles' volume fraction is seen from Figs. 12͑b͒ and 14 to be underestimated, while the large particles' volume fraction is overestimated. To decrease the relative importance of the attenuation at high frequencies, we solved a slightly different inverse problem in which both the attenuation and ␣ were divided by f 2 . However, only small improvements were found by modifying the attenuation data this way. The inverse-problem result shown in Fig. 14 was in fact obtained in this way.
Some insight into why the choice of f max drastically affects the results may be gained from Fig. 16, which shows the three-dimensional plots for the kernel K(a i ,a j ) for the same values of f max as considered in Fig. 12. We see that when f max ϭ10.4 MHz, corresponding to the circle in Fig.  12͑a͒, the kernel has a maximum for a i ϭa j ϭa max . The kernel for smaller particles is very small and, as a consequence, the inverse technique could determine the larger particle size volume fraction correctly but failed for smaller particles. In contrast to this the kernel for f max ϭ17.1 MHz, corresponding to the end of first peak, shows significant variations for a wide range of values of a i and a j , and this apparently leads to a much better inverse solution. Finally, the kernel for f max ϭ30.4 MHz, corresponding to the end of the third resonance peak, shows a less pronounced structure.
It is also instructive to examine the kernel and the results of inverse calculations for the problem of determining bubble-size distribution in bubbly liquids examined by Duraiswami. 2 The inverse procedure works very well for bubbly liquids as can be seen from Fig. 17͑a͒ which shows the input and computed bubble size distributions to be in excellent agreement. The kernel for this case has smooth variations over a wide range of bubble radii as seen in Figure  17͑b͒. The attenuation as a function of frequency is shown in Fig. 17͑c͒. The main reason for the success of the inverse technique for bubbly liquids seems to be that there is one resonance frequency for bubbles of each size. As long as the frequency range is broad enough to cover the resonance frequency of all the bubbles, it is possible to determine the size distribution.
The results presented so far were for polystyrene particles. We have also carried out inverse calculations for glass particles. As indicated earlier ͑cf. Fig. 6͒ there is no clear, sharp resonance frequency peak for glass particles. As a consequence, the inverse calculations for the glass particles did not show, in general, good agreement with the input size distribution.
The results presented so far show that the success of Tikhonov regularization to solve the inverse problem is limited. Although we have given plausible reasons for why the FIG. 14. As in Fig. 10͑b͒, but after cutting off the frequency range over which the attenuation was given between the first and second ͑attenuation͒ peak, indicated by a triangle in Fig. 10͑a͒.   FIG. 15. As Fig. 12, but with a broader size distribution. method works well for bubbles but not for all particles, it is possible that other techniques for solving the inverse problem may be more successful. For that reason we have attempted an alternative method 2,3,19 based on linear programing that we shall briefly describe here.
The constraint (a)у0 for all a was satisfied only a posteriori in the Tikhonov scheme. To ensure that the error is minimized while satisfying this constraint, we reformulate the original inverse problem as an optimization problem. The simplest scheme is to minimize the error

͑31͒
instead of the integral of the square of the quantity enclosed by two vertical bars at each frequency. Constraints on the solution are used a priori in optimization via linear programming; here we use that (a)у0. Imposing an upper bound on the total volume fraction ͑maximum packing͒ can also be incorporated but is not essential. After discretizing the frequency range by M and (a) in N discrete values we write ͚ jϭ1 N B i j ͑a j ͒Ϫ␣ tot ͑ f i ͒ϭu i Ϫv i , u i ,v i у0, iϭ1,2, . . . ,M.

͑32͒
Here, B i j is the discretized form of the integral operator in ͑31͒ and u i and v i are, as yet, unknown, non-negative variables. Now, it can be shown 19 that minimizing the absolute value of ͑32͒ is equivalent to minimizing with ͑32͒ as a constraint together with the constrains u i , v i у0 (iϭ1,...,M ) and (a i )у0 (iϭ1,...,N). Essential here is the notion that at the optimum u i v i ϭ0 for each i, which makes the solutions of the two minimization problems ͑31͒ and ͑33͒ identical. The above scheme was applied to a number of cases that were also examined using the Tikhonov method. It was found that, in general, the linear programing scheme produced inferior results. A typical example is shown in Fig. 18 where the Tikhonov method is seen to yield far better results for the size distributions. This technique also did not yield good inverse results for the case of glass particles.

V. CONCLUSION
A theory for the attenuation and wave speed of solidliquid suspensions at low particle volume fractions is described. The theory is shown to be in excellent agreement with the experimental data measured in our laboratory. Tikhonov regularization and linear programing techniques are employed to solve the inverse problem of determining the particle size distribution from the attenuation-frequency data. Although these techniques are successful in solving the inverse problem in several cases, we have also shown that the results are very sensitive to the choice of frequency range, the physical properties of the particles, and the nature of particle size distribution ͑unimodal, bimodal, etc.͒. Since the same techniques worked very well for bubbly liquids, we attribute the failure in solving the inverse problem satisfactorily to the complex resonance behavior of slurries. We conclude therefore that the prospects of using acoustic probes for on-line monitoring of particle size distribution of slurries are somewhat limited unless some additional information on the particle size distribution ͑e.g., unimodal͒ is available.
FIG. 16. The kernel K(a i ,a j ) for polystyrene particles when using for f max the value indicated in Fig. 12͑a͒ by a ᭺ ͑a͒, ϩ ͑b͒, and ᮀ ͑c͒. In both cases the attenuation was cut off at the same frequency, indicated by the square in Fig. 12͑a͒.