1 Introduction

The de facto standard non-invasive flow measurement techniques in research laboratories are based on optical principles (Tropea et al. 2007), the most prominent examples being laser-Doppler anemometry (LDA) and particle image velocimetry (PIV). While very successful for single-phase flows in e.g. wind tunnels or water flumes, the presence of a dispersed phase with a volume fraction as low as 0.5% can render these techniques useless (Deen et al. 2002; Poelma et al. 2006). The exact limiting volume fraction depends on the dispersed phase size, distribution and the dimensions of the flow domain (Linne et al. 2009). For small particles, droplets, and bubbles the limit is reached earlier, while it may be possible to obtain results ‘in between’ large bubbles with an overall large volume fraction (Mudde et al. 1997). The latter approach case may lead to biased statistics, however. For some studies, it is feasible to use refractive index matching, so that the system remains optically transparent (Wiederseiner et al. 2011). Unfortunately, only a limited number of combinations of solid and liquid materials can be used to achieve this matching, severely limiting the physical parameter space that can be probed. This means that many important two-phase flow problems, inspired by e.g. industrial or geophysical applications, are out of reach for the current optical measurement techniques.

In the last decade, a series of non-optical imaging modalities have been introduced for flow measurement, most of them based on medical imaging techniques. Examples include Magnetic Resonance Imaging (Elkins and Alley 2007; Ooij et al. 2011; Lakshmanan et al. 2016), X-ray imaging (Fouras et al. 2007; Heindel 2011) and ultrasound imaging (echography). Alternative strategies use tomographic reconstruction of electrical signals from sensors surrounding the flow, such as electrical capacitance tomography and electrical impedance tomography (ECT, EIT) (Yang and Peng 2002; Dyakowski et al. 2000; George et al. 2000). This review focusses on a particular subset of ultrasound techniques, namely correlation-based velocity estimation using echography image data; this is distinct from Doppler-based methods, such as e.g. Colour Doppler or Ultrasound Doppler Profiling. For clarity, these Doppler methods are also shortly described. For an extensive review of these techniques in a clinical setting, the reader is referred to the recent two-part review by Jensen et al. (2016a, b). Correlation-based techniques are known under various names in literature, including Ultrasound Imaging Velocimetry (UIV), speckle tracking velocimetry, and echo-PIV. They all refer to methods that apply PIV data processing techniques to data that has been obtained using echography (i.e. ultrasound imaging). Despite the alternative imaging modality, UIV retains the main strength of PIV: it produces non-invasive, instantaneous, two-dimensional, two-component velocity fields. An illustrative application example is shown in Fig. 1.

Fig. 1
figure 1

Reproduced from Gurung et al. (2016) ©IOP Publishing. Reproduced with permission. All rights reserved

An illustrative example of Ultrasound Imaging Velocimetry. (Left) a linear transducer placed at the wall of a pipe containing an opaque flow (here a strongly non-Newtonian drilling fluid). Reproduced from Poelma and Gurung (2016). (Right) the resulting vector field at the centerplane, superimposed on an instantaneous ultrasound image

Within the imaging modalities, ultrasound-based techniques represent a good compromise between imaging/measurement capabilities and practical considerations and constraints. For instance, it is significantly more affordable than MRI, which is also bound to a specific site (void of ferromagnetic materials) and which requires more complex postprocessing of data. Ultrasound-based techniques do also not require the stringent radiation regulations surrounding X-ray imaging. Furthermore, the latter is a shadowgraphy technique: the signal is integrated along the beam axis, complicating the analysis of three-dimensional flows. A tomographic approach is a possible solution (Fouras et al. 2007), but this is suitable for time-averaged flow patterns only. Optical coherence tomography (Huang et al. 1991) has a superior resolution, but is restricted to a small field-of-view. Finally, the resolution of ultrasound surpasses that of ECT/EIT, which generally uses a limited number of sensors. With this limited amount of information, reconstruction of complex phase distributions within the measurement domain becomes difficult, as the inverse problem at the heart of the method is severely ill-posed.

The accessibility, relative simplicity and low cost, combined with ever increasing possibilities, has resulted in a strong increase in the use of UIV in the last decade. Applications range from in vivo cardiovascular flow studies to more traditional fluid mechanics studies. In this manuscript, the basic concepts of ultrasound imaging (Sect. 2) and velocity measurement (Sect. 3) are reviewed. While UIV borrows heavily from PIV, there are some fundamental differences; these are described in Sect. 4. Examples of UIV applications are given in Sect. 5. Conclusions and some future opportunities are given in Sect. 6.

2 Ultrasound imaging

As this review is primarily intended for the fluid mechanics community, a brief historical overview (Sect. 2.1) and a basic introduction into ultrasound imaging (Sect. 2.2) will be given. Excellent in-depth discussions are found in e.g. Szabo’s monograph Szabo (2004) and various medical imaging textbooks, e.g. Suetens (2009).

2.1 Historical perspective

Ultrasound, defined as sound with a frequency above 20 kHz, has been used for diagnostics since the 1940s (Szabo 2004). While ultrasound was an established phenomenon in the nineteenth century, practical ways to harness it were enabled by the discovery of the piezo-electric effect by Jacques and Pierre Curie in 1880–1881 (Curie and Curie 1881). Research built on work in the field of SONAR, which emerged during the First and Second World Wars, and was used to locate submerged objects using echolocation. One of the earliest peacetime applications was an apparatus to detect flaws in metals—the ‘Supersonic Reflectoscope’ by Firestone (1946). The device made use of the property of ultrasound that it can penetrate non-transparent materials: first, a short pulse of ultrasound waves was created using a quartz-crystal and sent through a test sample. If a flaw was present in the sample, the variation in the material properties (in this case impedance) resulted in a reflection of a part of the sound waves. The intensity and time delay of these reflections, recorded using the same crystal and visualised on an oscilloscope, could then be used to localise and quantify the flaw. Around the same time, Dusik pioneered medical diagnostics using ultrasound (Dussik 1942). In his earliest experiments, he used a transmission approach (i.e. placing the receiving crystal beyond the field of interest) to study the human brain; this led to the first ever ultrasound image recorded.

At the core, modern echography still operates along the same lines as Firestone’s device. Modern devices use piezo-electric transducers to create and receive ultrasound. The received signals are naturally no longer directly displaced on an oscilloscope, but are digitised and processed further to form an image. Initially, only a single transducer was used. The 2D image created by Dussik was obtained by manually translating the single transducer (i.e. a point-by-point scan). Later, freehand scanning was introduced, which tracked the position and orientation of the transducer; by combining data from various viewing angles, this allowed reconstruction of stationary images in 2D and, later, in 3D (Barry et al. 1997). Mechanical scanners introduced faster imaging rates, so that (successive) 2D snapshots could be obtained. This was achieved by encasing a linear translation mechanism for the transducer element (Wild and Reid 1952). Alternatively, the transducer could be revolved or rocked around its axis, leading to sector imaging (with its characteristic wedge-shaped field of view). Developments in transducer materials, design and electronics led to phased-array probes (Bom et al. 1973; Ramm and Thurstone 1976). These phased-array transducer form the basis of modern equipment. They construct images by sequentially reading out elements, rather than mechanically translating an element. This has additional advantages, such as the possibility of beam steering and focussing. As a linear phased-array creates a 2D image, a 3D image can be constructed by sweeping the probe (Fenster et al. 2001). Alternatively, matrix transducers have been introduced (Smith et al. 1991) to obtain 3D data. Apart from such hardware developments, a better understanding of ultrasound wave physics and more sophisticated signal processing techniques have led to better image quality; a prime example is harmonic imaging (Duck 2002). In this method, the original frequency is filtered out and only non-linear responses are used for image reconstruction. This way, particular structures can be highlighted that would otherwise be obscured (e.g. contrast microbubbles; see later). Furthermore, the introduction of plane-wave imaging has increased the imaging rate significantly (Tanter and Fink 2014). In this technique all transducer elements are used simultaneously (rather than sequentially) and the image formation is done off-line using the recorded signals, analogous to digital holography (Tanter and Fink 2014; Schnars and Jueptner 2005). Apart from these technical breakthroughs, there have been many ultrasound innovations that are application-driven, such as the development of miniaturised transducers for intravascular ultrasound imaging (‘IVUS’, Roelandt et al. 1989) and specific transducers for specific organs and/or procedures (Szabo and Lewin 2013).

For this review, the focus is on conventional, linear phased-array transducers, as these are common in experimental fluid mechanics. In Sect. 6, the implications of some new developments on UIV will be discussed. A schematic representation of a linear phased-array transducer is shown in Fig. 2. Typically, these transducers have 64 up to 512 elements, which can emit and receive independently in specific patterns for a wide range of imaging modes. Each transducer has a certain inherent centre frequency (based on the resonance frequency of the particular piezo-electric element design), but can operate in a wider range (bandwidth), e.g. 5–10 MHz. Figure 2 also introduces the conventional coordinate system and terminology, which will be used and explained in subsequent sections.

Fig. 2
figure 2

Schematic representation of a linear phased-array transducer and coordinate system. ABW azimuth beam width

2.2 Ultrasound imaging basics

Figure 3 demonstrates the basic principles for the acquisition of an ultrasound image. Consider one transducer element of a phased-array. An alternating electrical signal is applied to this element, typically made of PZT (lead zirconate titanate) or PVDF (polyvinylidene fluoride) (Suetens 2009). The alternating signal leads to a deformation of the piezoelectric element and thus acts as a source of a pressure wave. Generally, relatively short pulses are created, with a frequency f. For medical applications, generally frequencies in the range of 1–10 MHz are used. High-frequency devices exist for specific applications, which operate in the 30–100 MHz range (Lockwood et al. 1996). For simplicity, it is here assumed that this pressure pulse is a small-amplitude, longitudinal compression/rarefaction wave travelling along the z (beam) axis, i.e. perpendicular to the transducer array.

Fig. 3
figure 3

Schematic representation of ultrasound imaging: a Overview of geometry, b raw RF signal, c log-compressed intensity signal, and d B-mode image. See text for details

The longitudinal wave travels with a velocity c, the speed of sound of the medium that carries the wave. If this wave encounters a medium with different properties, it may partially be transmitted and partially be reflected. The strength of the reflection is given by the difference in the specific acoustic impedance, as discussed in detail in the next section. In the schematic of Fig. 3a, the pulse first hits the surface between gel and the wall containing the flow. This creates a strong reflection, which is recorded by the same transducer element. Figure 3b shows the voltages (corresponding to pressure variations) from the transducer element as a function of time. At around 15 ms the strong echo from the wall can be observed. By using the speed of sound, the time axis of the collected data can be converted into a spatial axis, as shown in Fig. 3c. This type of data is called A-mode, with ‘A’ referring to amplitude. Several additional steps were performed going from (b) to (c): a Hilbert transform is used to remove the fast modulation on top of the signal (also referred to as envelope detection), see also the small inset in subfigure (b). Furthermore, the data is log-compressed, to reduce the difference in intensity of specular reflections (e.g. from flat surfaces) and scattered signals (e.g. blood, tracer particles).

Once the data from this particular transducer element has been collected, the process can be repeated using the subsequent element. The result is a 2D ultrasound image, as shown in Fig. 3d. The dashed vertical line represents the data shown in subfigure (c); each vertical line is also referred to as a scan line. These 2D results are referred to as B-mode images, with b referring to ‘brightness’. Alternatively, the same transducer element can be read out again and again. This creates M-mode data, where the vertical axis represents the z position (distance from the transducer) and the horizontal axis represents time. This mode is often used to study fast dynamic effects in the body, such as e.g. the motion of heart valves.

Depending on the conditions of the experiments, time gain compensation can be used. This compensates for attenuation of the signal as it travels deeper and deeper into the imaged region. Generally, the details of the correction function are manually specified by the operator, in order to equalise the intensity across the image.

As a final step, the acquired image may need to be corrected for the shape of the transducer. A typical example of this so-called scan conversion is shown in Fig. 4. Here the raw data B-mode data (i.e. an intensity line per transducer element) is transformed to the corresponding wedge-shaped image that represents the true physical field-of-view for the curved transducer that was used. Naturally, this process requires some sort of interpolation process that leads to an non-uniform spatial resolution and to distorted particle images.

Fig. 4
figure 4

Reproduced from Kim et al. (2004a). ©Springer-Verlag 2004. Used with permission

Raw B-mode data of microbubbles in water (left) and the corresponding scan-converted (right) image

Implicitly, it has been assumed that the emission and receiving was done sequentially using individual elements. In practice phased-arrays are used in more elaborate ways, generally using groups of elements to record a single A-mode line. A prime example is focussing, as is also demonstrated in Fig. 2: here five elements are used, triggered with an electronic delay so that the wave fronts combine to a focal region. The number of elements that is used determines the aperture. Naturally, this method of focussing provides more flexibility than a fixed, physical acoustic lens. Note that such a lens is still used to focus the beam in the elevational direction (as opposed to the azimuthal focussing done by the phased-array). In Fig. 2 the focal location in the elevational and azimuthal plane are identical, but this is not required. Multiple focal points at various depths can often be selected, but this comes at the cost of imaging rate, as they are acquired sequentially. Other uses of a phased-array are beam steering (i.e. tilting the beam axis), and a variable line density (with the option to create more lines than transducer elements by means of interpolation).

While the focus of this review is on velocimetry, it should be pointed out that ultrasound imaging by itself can also be a useful tool in fluid mechanics (beyond the obvious cardiovascular applications). An illustrative example is shown in Fig. 5, where the deformation of a compliant tube due to a traveling pressure wave is visualised using echography. This example, taken from Hickerson and Gharib (2006), demonstrates that imaging is possible without the need for refractive index matching, which can be difficult to achieve if also mechanical properties (e.g. Young’s modulus) of the wall need to be matched.

Fig. 5
figure 5

Image reproduced from Hickerson and Gharib (2006). ©Cambridge University Press. Used with permission

A pressure wave traveling through a compliant tube (from left to right) is visualised using two consecutive ultrasound images. The position of the field-of-view and transducer location is shown schematically

2.3 Reflection, refraction and scattering

Ultrasound is an archetypical wave phenomenon. The (small) pressure fluctuations satisfy the wave equation and thus typical phenomena like wave reflection, refraction and scattering can be observed. While in optics the refractive index is the key parameter, for (ultra)sound this parameter is the specific acoustic impedance Z, the product of the speed of sound and the density,Footnote 1 \(Z\equiv \rho c\).

Consider a small-amplitude planar wave, with amplitude (or intensity) \(I_i\), hitting a surface with an incidence angle \(\theta _i\), see Fig. 6. The surface separates the two regions with impedance \(Z_1\) and \(Z_2\) and speed of sound \(c_1\) and \(c_2\). The angles of the reflected (\(\theta _{\text{r}}\)) and refracted (‘transmitted’; \(\theta _{\text{t}}\)) waves follow from Snell’s law:

Fig. 6
figure 6

Snell’s law

$$\begin{aligned} \frac{\sin\ \theta _i}{c_1} = \frac{\sin \ \theta _t}{c_2} = \frac{\sin \ \theta _r}{c_1} \end{aligned}$$
(1)

The ratio of the amplitudes of the refracted and reflected wave, the reflection coefficient R, can be calculated from the Fresnel equations (Suetens 2009):

$$\begin{aligned} R \equiv \frac{I_r}{I_i} = \frac{Z_2\ \cos \ \theta _i - Z_1\ \cos \ \theta _t}{Z_2\ \cos \ \theta _i + Z_1\ \cos \ \theta _t} \end{aligned}$$
(2)

For a planar wave that hits a surface with \(\theta _i=0\), as was the case in the schematic representation in Fig. 3, this reduces to:

$$\begin{aligned} R = \left( \frac{Z_2-Z_1}{Z_2+Z_1} \right) \end{aligned}$$
(3)

Note that the reflection coefficient is sometimes defined with a squared right-hand-side. In this case the coefficient described the reflected power, rather than amplitude. The transmission coefficient can be defined as \(T=1-R\). When both media have matched impedances (\(Z_1 = Z_2\)), the reflection coefficient is zero and the wave is transmitted completely (\(T=1\)). A negative reflection coefficient implies inversion of the wave shape.

Table 1 provides a list of materials and their acoustic properties, with a focus on materials used to construct fluid mechanics facilities. Note that this data is intended to be indicative; different conditions and the exact material properties (e.g. molecular weight of polymers) may lead to slightly different values.

Table 1 Acoustic properties of common materials

Most medical ultrasound textbooks provide extensive tables for various tissue types found in the body (liver, kidney, blood, etc.). A significant observation is that the speed of sound within the body is confined to a fairly narrow range (1510–1584 m/s), apart from bone (3198 m/s). As the speed of sound varies less than 5%, the conversion of time delay to a spatial position (i.e. the step from Fig. 3b, c) can be achieved by a simple multiplication. In other words, there is a single calibration constant in the z direction to convert the data to physical units. For systems containing more variation, such as strongly stratified flows or facilities with thick walls, this would lead to distorted images and a more elaborate calibration is required to account for the local speed of sound. Note also that temperature has a strong effect on the speed of sound: for water, for instance, an increase of 10 °C leads to a difference of a few percent. It is thus wise to perform an in situ check (e.g. by comparing the observed diameter of a tube with the nominal value) to ensure accurate scaling of the image.

From the data in Table 1 and Eq. 3, it can be seen that the common materials that are used for the construction of test sections (based on experience with optical techniques) are far from optimal for ultrasound imaging. Generally, the transducer is covered with a layer of coupling gel to ensure a good contact with the object under investigation. The transition from this gel to e.g. glass will lead to a reflection of nearly 40% of the amplitude; subsequent transitions (glass to water) will further reduce the signal strength. For comparison: the transition from air to glass will only reflect 4% of an incoming ray of light. Several recent studies have used materials and test sections that are tuned for ultrasound imaging, rather than optical techniques (Gurung et al. 2016; Gurung and Poelma 2016). There is no need for optical transparency of the material, but rather a low impedance in combination with appropriate mechanical strength. Neoprene and polyethylene are a good choice in combination with water if ultrasound image quality is a concern.

In optics, thick-film antireflective coatings reduce the amount of reflected light by using a layer with a refractive index that is in between the values of air and the lens material. A similar approach is used in the design of ultrasound transducers, where ‘matching layers’ ensure that the acoustic output from the piezoelectric materials is maximised. This could also be used in the design of experimental facilities to overcome large jumps in impedance (e.g. from gel to glass or metal). By placing an additional layer with an intermediate acoustic impedance, the percentage of transmission can be increased.

From Table 1, it can also be concluded that the transition from a solid or liquid to air (or vapour) leads to a near-complete specular reflection of the signal. This implies that imaging through air cavities or large bubbles is not possible. Nevertheless, ultrasound techniques can still be used to reconstruct relatively simple shapes. For instance, Felici et al. (2013) used the strong echo from the water / vapour transition to characterise the shape of a cavitation bubble on a hydrofoil. The strong reflection occurring at the transition from liquid to gaseous media also explains the need of placing contact gel on the transducer surface. Similarly, it explains why imaging of the heart is often done endoscopically via the esophagus (‘transesophageal echocardiogram’), to avoid having to image through the lung cavities. Note that insoluble liquids (Castor oil, turpentine, etc.) generally have a slightly lower impedance than water. While this may lead to interesting applications in the study of emulsions, so far no studies have reported imaging in these systems.

The pressure pulse propagation is not a lossless process, even in homogeneous media. Attenuation of the signal can occur due to reflection, refraction, spreading loss (for spherical waves), scattering, and absorption Szabo (2004). The latter is due to partial conversion of the pressure pulse to heat, due to viscosity in fluids and viscoelasticity in solids. For Newtonian fluids, the attenuation coefficient can be predicted theoretically using Stokes’ law of sound attenuation. Attenuation by solids is often modelled as an exponential decay, using an empirical absorption coefficient \(\alpha\). Absorption is strongly dependent on the wavenumber, which is reflected in the units for \(\alpha\): for instance, water has an absorption coefficient \(\alpha = 2.17\times 10^{-3}\, {\text{dB}}/({\text{MHz}^2}\, {\text{cm}})\), while that for blood is \(\alpha = 0.15\,{\text{dB}}/({\text{MHz}^{1.21}}\,{\text{cm}})\) (Szabo 2004). The non-linear nature of absorption results in a change of the pulse shape, not only its amplitude. Generally speaking, the absorption is inversely proportional to the frequency for various tissue types (note the exponent in the units of the values of \(\alpha\) reported above). As the resolution is determined by the frequency, this leads to a trade-off between imaging depth and image resolution. As mentioned in Sect. 2.2, the signal attenuation is generally countered by applying time gain compensation, an electronic amplification of the received signals to obtain a constant signal strength for all depths. The effectiveness of this procedure depends on the dynamic range of the hardware.

As stated, scattering is a major contributor to the attenuation of the signal. This process is different than the specular reflections that are created by e.g. walls, i.e. objects that are much larger than the wavelength of the pulse. In contrast, smaller objects scatter omni-directionally, so a much weaker signal is received by the transducer. The exact scattering properties are determined by the size of the scattering object with respect to the wavelength of the sound source. Isolated objects much smaller than the wavelength act as point scatterers (diffusive or Rayleigh scattering). The scattered intensity is strongly dependent on the wavelength (\(I \sim f^4\)), unlike specular reflections (Szabo 2004). Note that for ultrasound imaging, the amount of backscattered waves (\(\theta = \pi\)) is the most relevant, unlike e.g. conventional PIV applications (were typically \(\theta = \pi /2\) is the most relevant). For objects that are of the same size as the wavelength, the diffractively-scattered signal may reach the transducer with slightly different path-lengths (e.g. from the ‘top’ and ‘side’ of an object). This introduces interference and thus a speckle-like signal is received. The same holds for densely concentrated point scatterers. Heterogeneous materials (such as tissue, which has continuous, minor variations of its acoustic properties) also exhibits this behaviour. While the theoretical description is the same as that for light interacting with particles, the wavelengths are much larger (compare the typical wavelength of 500 nm for light with \(1500/10\times 10^6 = 0.15\, \hbox {mm}\) for sound waves of 10 MHz in water). A typical flow tracer of 10 micrometer will be small compared to this wavelength and thus act as a Rayleigh scatterer. In contrast, typical ‘dispersed phase particles’, such as sand, will have a diameter comparable to the wavelength of ultrasound. In this case, the Rayleigh approximation no longer holds and the full Lorenz–Mie equation has to be solved to predict the theoretical scattering behaviour. The interested reader is referred to e.g. Szabo [for a description in terms of sound (Szabo 2004, Chapter 8)] or classic texts on optics (Born and Wolf 1999).

It is tempting to overcome the attenuation by simply increasing the power of the sound waves. However, there are three main practical issues that have to be taken into account: heat, cavitation and acoustic streaming. A considerable part of the electrical power is converted to heat in the transducer (Saunders et al. 2004). This so-called self-heating may be harmful for patients, but is also relevant for non-medical applications, as high temperatures may damage the transducer itself. The second limitation arises from the nature of the ultrasound waves: they are compression/rarefaction waves. When the (acoustic) power is increased, the increased amplitude will lead to lower values of the pressure in the rarefied parts of the wave. If the pressure drops below the vapour pressure, cavitation may occur. To avoid this risk in clinical applications, the mechanical index was introduced (Abbott 1999). This index is defined as the ratio of the peak negative pressure and the square root of the frequency (\(P_{\min}/\sqrt{f}\)); it cannot exceed 1.9 for clinical applications.

Finally, intense acoustic pressure waves may induce a mean motion of the fluid. While the underlying physics are well documented (Lighthill 1978; Duck et al. 1998; Sarvazyan et al. 2010), many textbooks give over-simplified and often contradictory explanations. There are three mechanisms that can lead to a mean motion of the wave-carrying medium: (1) finite-amplitude effects of the longitudinal wave give rise to a mean transport of fluid elements (cf. Stokes drift). This is generally a small effect Lighthill (1978). (2) If there is a change in impedance in the medium (density and/or speed of sound) an acoustic radiation force is created due to the partial reflection of the pressure wave. This will create a gradient in the energy density, which will push, for instance, suspended particles or microbubbles in the direction of the wave propagation (Dayton et al. 2002). (3) The same force also occurs in a homogeneous medium, if there is absorption. Attenuation of the acoustic energy flux in the fluid (by e.g. viscosity) introduces a Reynolds stress in the direction of the wave propagation. This Reynolds stress pushes the fluid away from the transducer. Note that attenuation in a medium can also be due to the presence of many small scatterers, so that the second and third mechanism may appear similar; however, in the second mechanism the force acts on the particles, while in the third it acts on the fluid. For dense suspensions of small particles this distinction may be academic, as the observed effect on the flow is the same. The resulting phenomenon caused by the acoustic radiation force due to absorption, i.e. a mean motion of the medium in the direction of sound propagation, is generally referred to as ‘acoustic streaming’ (Sarvazyan et al. 2010). An example of this effect can be seen in Fig. 7, where the same flow is measured using a low and high power acoustic field. Supplemental movie S1 also shows an example, where a steady flow is created solely by acoustic radiation forces. Various analytical models and approximations are available in dedicated texts [e.g. Lighthill (1978)] for either an unbounded medium or in the presence of a solid boundary. For practical situations, a more pragmatic approach is to obtain velocity fields using different power levels and compare the mean flow pattern (as done in Fig. 7).

Fig. 7
figure 7

An example of the effects of acoustic streaming. A laminar flow (from left to right) is observed using ultrasound with low and high power. The background image is generated by superimposing a series of images to visualise the flow using streaks. The colour-coding of the vectors of the PIV result highlights the secondary flow (in the beam axis direction); here the secondary flow for the high-power case is approximately 10% of the centerline velocity

2.4 Field-of-view, resolution and acquisition rate

The field-of-view that can be obtained depends on the transducer geometry and acquisition settings. The geometry determines the maximal width (\(L_x\), see Fig. 2; also referred to as ‘footprint’), assuming that all elements are used (‘100% sector’). For linear transducers, the field-of-view is generally rectangular, although beam-steering techniques can distort this to a parallelogram or wedge shape. Curved (‘convex’) transducers and mechanically rotating or rocking transducers have a wedge-shaped field-of-view after scan conversion, see e.g. Fig. 4.

The depth of the image (\(L_z\)) is specified during the acquisition. This depth corresponds to the maximum allowable time of flight (\(2\times L_z/\hbox {c}\)), i.e. until how long after emitting a pulse should echoes be recorded. The selected depth directly affects the acquisition rate (see later), but can also be a convenient way to block spurious echoes from beyond the region of interest.

The thickness of the measurement plane (\(L_y\), comparable to the light-sheet thickness in conventional PIV) is often referred to as the elevational beam width or slice thickness (Goldstein 1988). It is determined by the physical dimension of the transducers and the acoustic lens that is affixed to the array; recall that this lens only focusses in the elevational direction. As a result, the elevational beam width is not constant, but has a set minimum at some distance from the transducer. Typically, the minimum thickness is a few millimeters. Several methods have been introduced to quantify the effective measurement volume: for instance, the volume can be sampled point-by-point using a needle hydrophone to determine the local acoustic pressure (Poelma et al. 2012). Goldstein (1988) introduced a test object with a reflecting surface at 45°; the ‘height’ of the scattered band in the image can then be used to find the out-of-plane thickness. Alternatively, a small scatterer can be translated in the out-of-plane direction, while recording its image intensity (Gurung and Poelma 2016). For detailed investigations it is possible to apply Schlieren techniques to visualise the ultrasound wave field (Neumann and Ermert 2006).

The resolution of an ultrasound image is generally different in the scan (x, lateral) and beam (z, axial) directions. In each case, the resolution represents the smallest separation between two scatterers that can still be distinguished. To understand the axial resolution, consider the two points \(P_1\) and \(P_2\) in Fig. 8a, which are both on the same beam axis and separated by a small distance, \(\delta \hbox {z}\). A pulse is emitted and as soon as this reaches the scatterer at \(P_1\) an echo is created; for simplicity we assume this is an exact copy of the original pulse. The original pulse will reach the second scatterer (\(P_2\)) slightly later. The two echoes will be received with a temporal separation equal to \(2\delta z/c\). This implies that the pulse duration should be smaller than this temporal separation, as the echoes would overlap otherwise.Footnote 2 The pulse duration is given by the wave period (1/f) and number of cycles (n); alternatively, this can be expressed as a spatial length of \(\hbox {n}\lambda\). While a lower number of cycles may lead to a higher resolution, it comes at the price of a higher bandwidth. Transducers have a limited bandwidth, typically 0.2–0.5 times the centre frequency. A very short pulse (e.g. one sinusoidal cycle) cannot be captured efficiently. On the other hand, long pulses have a narrower bandwidth, but lead to a lower resolution. Generally, 2–4 cycles are used in pulsed ultrasound as a compromise between bandwidth restrictions and resolution. Using a pulse length of \(n=2\), the axial resolution is found to be equal to the wavelength \(\lambda\): \(\delta \hbox {z} = \lambda = c/f\). For an ultrasound frequency of 5 MHz in water, this gives an image resolution in the beam direction of 0.3 mm. Higher frequencies will give higher resolutions, but at the cost of more attenuation, as discussed in the previous section. Note that RF data is usually sampled at a higher frequency than the original ultrasound frequency (e.g. sampled at 40 MHz for a transmitted frequency of 10 MHz). This sampling rate, converted to its spatial equivalent, should not be confused with the axial resolution.

Fig. 8
figure 8

Schematic illustration of axial and lateral resolutions. a A short pulse encounters two scatterers, separated by a distance \(\delta \hbox {z}\) along the beam axis; b the wave field of a single element; c the wave field of 5 elements firing with specific temporal delay to create a focal point. For clarity continuous waves are used in (b) and (c), rather than pulses

The resolution in the scan (lateral) direction is less straightforward. Several factors play a role: the size of a single transducer element, the pitch (the distance between the centres of two adjacent elements), wavelength/frequency and the beamforming method (as discussed below). Each transducer element emits a beam with an initial size close to the dimensions of the element. This beam is focussed in the elevational direction and diverges in the azimuthal plane, see also Fig. 8b. As stated before, the use of a group of elements in a phased-array allows the creating of a focal region (in the azimuthal plane) by means of beamforming. This is a technique where signals from various elements are combined (with appropriate delay) to create a directed and/or focussed wave front. This is shown schematically in Fig. 8c, where 5 identical sources—the same as in panel (b)—are used. The first and fifth element fire first, followed by the second and fourth and finally the central element. The interference of these wavefields creates a focal region where the lateral resolution is optimal. For an effective aperture size D, focal length F and wavelength \(\lambda\) the azimuthal beam width (ABW) at the focal point is \(\delta x = F\lambda / D\). If an image is constructed by merging data using focal points at various depths, the resolution will change with depth if the aperture is not adjusted accordingly.

The phased-array approach using multiple elements can also be used to obtain additional scan lines, i.e. images can be created with more scan lines than transducer elements. Keep in mind that the lateral resolution will still be determined by the previously mentioned parameters (effective aperture, focal length) and the image acquisition rate will decrease with an increasing scan line density. Typical resolutions in the lateral direction that can be achieved are in the range of 0.3–0.5 mm. The ‘out-of-plane’ resolution, \(\delta y\), can be estimated from the local elevational beam width. A convenient way to describe the effect of the finite resolution of a system is to report the observed image of a point scatter as an ellipsoid, with each of the axes representing the relevant resolution (axial, lateral, azimuthal) (Szabo 2004). The volume of this ellipsoid represents one ‘resolution cell’, which will be a useful concept to classify the nature of the images. A more formal way to define the resolution is to determine the point-spread-function (PSF) of a small scatterer (Contreras Ortiz et al. 2012); for practical purposes the resolution cell and PSF are generally equally useful.

Apart from the dimensions and resolution of the images, it is also essential to consider the image acquisition rate (‘frame rate’, FR). This is of particular importance for the correlation-based velocimetry methods described in the next section. In most modern devices the acquisition rate is no longer limited by electronic bandwidth, but purely due to time-of-flight of the pulses. A pulse needs \(\tau = 2L_z/c\) to reach the bottom of the field-of-view and return to the transducer. In water, and assuming a 4-cm imaging depth, this requires approximately \(54\,\mu \hbox {s}\). While this may not seem excessive, it is important to reiterate that in ultrasound imaging the image is build up from sequential scan lines. Constructing an image with 128 scan lines would require \(128\tau =6.9\,\hbox {ms}\). Assuming that there is no overhead in the system, this yields a frame rate limit of 145 fps. For imaging purposes this may be fast enough for many events. However, for velocimetry it sets a limit on the maximum velocities that can be obtained: The frame rate directly dictates the delay between successive frames (\(\Delta T = 1/{\rm FR}\)). \(\Delta T\) is a key parameter to optimise PIV measurements, in particular to ensure a moderate tracer displacement within the entire field of view. As the magnification is fixed in ultrasound imaging and most transducer designs have comparable geometries (in particular pitch), this limitation results in a practical limit of 0.5–0.8 m/s reported in many studies; the maximum reported velocity in many validation experiments is generally inversely proportional to the scale of the experiment, i.e. the imaging depth and the diameter of the tube (Liu et al. 2008; Poelma et al. 2011, 2012; Westerdale et al. 2011; Walker et al. 2014).

A solution to increase the acquisition rate is to reduce the imaging depth and/or number of scan lines. The latter can be done by e.g. only reading out a central subgroup of the transducer array. As an extreme example, Beulen et al. (2010) only read out 14 elements to obtain a frame rate of \(730\,\hbox {s}^{-1}\). Naturally, this reduces the lateral field-of-view, but in many cases even a single velocity profile (in the beam direction) will already provide sufficient information about the flow. Alternatively, one can improve the acquisition rate by skipping elements (i.e. a scan line density below 100%), but this will impact the resolution.

Several more elaborate approaches have been introduced to improve the acquisition rate. In ‘parallel beamforming’, several subgroups create separate beams simultaneously (Shattuck et al. 1984). The main idea is that these beams are not overlapping and thus several scan lines can be obtained simultaneously. ‘Synthetic aperture imaging’ uses only a single element to insonate the entire field-of-view, while still receiving with all elements. This comes at the cost of a much lower resolution. More recently, ‘plane wave imaging’ made it feasible to obtain images at kHz rates (Tanter and Fink 2014). All elements emit at the same time, creating a plane wave in the beam direction. Beamforming is done in the receiving stage (i.e. as a data post-processing step), again this comes at the price of a lower resolution and signal-to-noise ratio and requires specialised hardware. The previous methods primarily originate in the need for fast medical imaging, e.g. to perform elastography—which requires the measurement of fast shear waves in tissue. Recently, a method was introduced specifically for correlation-based velocimetry: interleaved imaging (Poelma and Fraser 2013). In this method, an image is still recorded by reading all elements sequentially. However, one image frame is recorded during ‘odd’ time-steps, while the ‘even’ time-steps are used to record a second image, with a particular off-set. A typical read-out sequence is shown in Table 2 to illustrate this process and an example is shown in Fig. 9. In this method, the total time for one image pair is the same as in a conventional method; however, the delay between the first and second image of a pair is flexible. Special attention is given here to this technique because it is a prime example of a data acquisition approach that is optimised for flow measurement, rather than for imaging. It is to be expected that many other new acquisition methods will appear in the future, inspired by progress in the medical ultrasound community.

Table 2 Normal and interleaved imaging scan line order example
Fig. 9
figure 9

Image reproduced from Poelma and Fraser (2013). ©IOP Publishing. Reproduced with permission. All rights reserved

Interleaved imaging example. The left-most images shows the scan lines as a function time; a close-up is shown in the middle to show the interleaved image, with an off-set of m timesteps. The right-hand image shows the de-interleaved images. The yellow dot is at a fixed location to show the tracer displacement

2.5 Hardware

A wide range of ultrasound imaging systems is available (Powers and Kremkau 2011). They range from portable systems, which fit in a suitcase, to sophisticated research-oriented systems. Recently, systems have been introduced that connect a scan head to a mobile phone, making them ultraportable. The associated cost for these systems ranges from a few thousands to hundreds of thousands euros for top-of-the-line systems. A lively market in refurbished hardware exists, which may be suitable for fluid mechanics laboratories, who do not need to adhere to the strict quality assurance regulations or maintenance cycles of hospitals.

State-of-the-art systems are available from most major healthcare technology companies (e.g. Siemens, Philips, GE).Footnote 3 These are fairly flexible and can be used with various transducer types and operate at a range of frequencies. Apart from imaging, additional capabilities are generally available for e.g. velocimetry and advanced diagnostic techniques. It should be noted that these devices are developed for robust, clinical use. This means that they are generally closed systems using proprietary algorithms and software. It is often difficult (or even impossible) to access raw RF data. It may not be evident what kind of processing (compression, filtering) occurs when image sequences are exported. This makes them less than optimal for flow measurement. Additionally, patient safety regulations may set limits to certain parameters (e.g. mechanical index). Some of these companies provide access through hardware and/or software add-ons that partially open up the system.

Several companies (also) produce systems that are more research oriented, such as e.g. ESAOTE, TELEMED and BK Ultrasound (formerly Ultrasonix). They provide more flexibility, often with community-driven software tools, which allows easier experimentation with new acquisition protocols. These systems represent a good compromise between user-friendliness and measurement capabilities, as well as access to raw RF data.

Systems that cater primarily to the research community are also available. Examples include ULA-OP (Tortoli et al. 2009), a commercially-available open research platform developed at the University of Florence (Italy), the Vantage platform (Verasonics; Kirkland, WA, USA) and the ‘OPEN’ system by Lecoeur Electronique (Chuelles, France). Finally, a series of specialistic, experimental devices have been developed, such as e.g. the RASMUS and SARUS systems at the Technical University of Denmark (Jensen et al. 2005, 2013) and the system developed specifically for echo-PIV by Liu et al. (2008). Such devices are generally developed in close collaboration between electronic engineering and medical imaging departments and have a high barrier of entry (costs, ease of use, etc.).

3 Velocimetry

Flow measurement using ultrasound has evolved concurrently with the imaging techniques discussed in the previous section. Ultrasound flow measurement can be classified in three groups: (1) time-of-flight methods, (2) Doppler methods and (3) imaging/correlation-based methods. The latter is the main topic of this review, but for clarity the other two are briefly discussed.

3.1 Time-of-flight methods

A prime example using the time-of-flight method is the ultrasonic flow meter, which is commonplace in process technology applications (Lynnworth and Liu 2006). When a pulse (transmitted at an angle with respect to the main flow direction) travels with the flow, the time-of-flight is shorter than in a non-moving fluid, where the time-of-flight would be dictated solely by the speed of sound of the medium. By using one emitter and two receivers (placed upstream and downstream on the opposite side of the flow), the mean flow velocity can be obtained without requiring knowledge of the speed of sound. The flow does not need to contain acoustic scatterers and the emitters/receivers can be placed non-invasively at the exterior of e.g. a pipe that is under investigation. A major drawback is that it measures the mean velocity along the pulse path and assumptions about the velocity profile need to be made to obtain flow rates. By using multiple transmitters, it is possible to reconstruct complex flow patterns using a tomographic approach. Since the early 1980s, this approach has been used to determine ocean currents (Behringer et al. 1982; Dushaw 2014). Manneville et al. (1999a, b) demonstrated that the technique can also be used in laboratory experiments. By using two arrays with 64 transducer elements on opposing sides of the region of interest, they studied the dynamics of a vortex. As the change in the time-of-flights is generally relatively small in a laboratory setting, they introduced a method to amplify the differences. This is done by repeatedly re-sending the time-inverted collected signals, so that the time-of-flight differences amplify with each pass (Manneville et al. 1999a). With seven iterations this approach improved the resolvable velocity considerably, reducing it from 2 down to 3 mm/s in their case.

3.2 Doppler methods

Doppler-based methods measure the change in an ultrasound pulse as it is reflected by an object that has a velocity component in the direction of the beam. This implies the fluid must contain scattering material, unlike the time-of-flight method. Note that the actual Doppler shift itself, i.e. the change in frequency, is generally not measured directly. Instead, the phase shift is determined between subsequent pulses originating from the same scatterer (Evans et al. 2011). This means that multiple pulses are required from one sample position. By selecting different delay times (‘time-gating’), the velocity at different depths can be obtained. This is illustrated in Fig. 10, based on the work by Takeda (1986). This paper introduced ultrasound Doppler techniques outside cardiovascular applications (Takeda 1986). The technique is generally referred to as ‘ultrasound velocity profiling’ (UVPFootnote 4) for implementations outside the medical field. Applications range from hydrological studies Gartner (2004) to in-line rheology in food production (Wiklund et al. 2007). A major advantage is that all electronics can be contained in a single, robust probe or pod and it does not need calibration. A major drawback is that only the velocity component along the beam axis can be obtained.

Fig. 10
figure 10

Illustration of the principle of Ultrasound Velocity Profiling: (a) configuration of ultrasound beam and flow, (b) received ultrasound signal and (c) reconstructed velocity profile. Redrawn based on Takeda (1986)

In the medical field, blood flow has been measured using the Doppler shift since the late 1950s, first using continuous signals and later using pulsed ultrasound (Sigel 1998). An important step forward was to combine (B-mode) images with the velocity data, which was first reported by Barber et al. (1974). This method is now referred to as ‘Colour Doppler’ or ‘Color Flow Imaging’: local blood velocity in a small region is colour codedFootnote 5 and superimposed on a grayscale image showing the anatomical information using B-mode images. As stated before, only the velocity component in the beam direction can be obtained. As blood vessels are generally parallel to the skin, this means that the dominant flow component cannot be measured easily. To overcome this, the beam direction is usually tilted by means of beam steering; this leads to a parallelogram-shaped Doppler region in the composite image. By assuming that the flow is predominantly unidirectional and aligned with the walls it is possible to estimate the streamwise velocities.

A drawback of Doppler velocimetry is aliasing, where local velocities lead to a phase shift that exceeds the detectable range set by the pulse repetition frequency (PRF). This PRF determines the time between successive pulses that are used to determine the phase shift (Evans et al. 2011) and needs to be optimised so that aliasing is avoided, while low velocities can still be resolved—analogous to the role of \(\Delta T\) in PIV/UIV. The drawbacks of angle-sensitivity and aliasing are largely overcome when using ‘Power Doppler’ (Rubin et al. 1994). This mode uses the power, i.e. the integral of the Doppler signal, rather than the mean frequency shift of the signal. A drawback is that it only shows the magnitude of the local velocity. Nevertheless, the excellent temporal resolution makes Power Doppler suitable for so-called functional imaging: it can visualise changes in blood flow patterns, e.g. in the brain as a result of external stimuli (Macé et al. 2011). Other techniques, such as MRI, are generally too slow to detect these changes with a sufficient spatial resolution. Recent work has combined ultrafast plane-wave and tomographic techniques, leading to ‘4D’ measurements in a rat brain, with a resolution of 100 μm × 100 μm × 100 μm and 10 ms (Demené et al. 2016).

Various ‘Vector Doppler’ techniques have been introduced to obtain both in-plane velocity components, rather than just the component along the beam axis. An example is the cross-beam approach, which uses two transducers placed at opposing angles (Dunmire et al. 2000). A similar multiple-angle approach can also be implemented with a single array using beam-steering (Fadnes et al. 2015). Another technique that was introduced is the transverse oscillation method (Jensen and Munk 1998; Jensen et al. 2011). In this method the ultrasound pulse (as shown in e.g. Fig. 8a) is modulated in the lateral direction. This modulation, on top of the original pulse shape, creates a pulse with a two-dimensional intensity pattern. The velocity component in the beam direction is obtained using conventional Doppler methods. The lateral component is obtained by cross-correlating the signals from adjacent transducer elements: a lateral displacement of the scatterer component will shift the modulated pattern. An application example using this technique is shown in Fig. 11. This approach gives instantaneous, two-component velocity fields in a manner very similar to the methods discussed in the subsequent section.

Fig. 11
figure 11

Image adapted from Udesen et al. (2007). ©2007 World Federation for Ultrasound in Medicine & Biology. Used with permission

‘Vector velocity imaging’; bifurcation of the carotid artery at the time of the peak systole

3.3 Correlation methods

Similar to optical PIV, which originated from speckle velocimetry in solid mechanics, UIV flow measurements where preceded by applications that estimated motion and deformation based on tracking of speckle patterns in tissue (Meunier and Bertrand 1995). Before addressing how motion—including flows—can be quantified, it is useful to address the nature of the images under consideration.

In most ultrasound images speckle is apparent, see e.g. the tissue outside the artery in the bottom left of Fig. 11. Tissue contains many microscopic (i.e. subresolution) scatterers. However, the observed overall pattern appears much more coarse-grained and is the product of the imaging process. The same holds for flows containing many small particles, as can be seen in the right-hand-side of Fig. 1; this fluid contains a 25% volume load of micrometer-sized particles, so that many scatterers are present within the volume representing the resolution. The intensity variations observed in the figure are due to speckle from interference, and do not represent individual particle images. In contrast, in Figs. 4 and 9 individual scatterers can be distinguished. Note that both experiments use micrometer-sized bubbles as tracers, yet their images are distinctly larger (order mm). For experiments with (very) low concentrations, tracking of individual scatterers is an option (Ackermann and Schmitz 2016). Using ultrafast imaging (>500 fps), it is even possible to localise the signal of individual microbubbles with sub-wavelength resolution (Errico et al. 2015).

In conventional PIV, images are classified based on their source density, which represents the number of scatterers per resolution unit (Adrian and Westerweel 2010). An equivalent expression can be formulated for ultrasound images:

$$\begin{aligned} N_{\text{S}} = \frac{C}{\delta _x \delta _y \delta _z} \sim \frac{C}{\lambda ^2} \end{aligned}$$
(4)

In this expression, C is the concentration of scatterers. The denominator contains the three resolution axes spanning a resolution cell as defined in Sect. 2.4. Substituting the expressions \(\delta _x = F\lambda /D\) and \(\delta _y = n \lambda\) (discussed in the same section), we find that the source density is proportional to the inverse of the wavelength squared. Using typical values (0.5 mm × 0.5 mm × 4 mm), we obtain a simple heuristic: \(N_{\rm S} \approx C'\), with \(C^{\prime}\) the concentration of scatterers per mm3. Images obtained with a low source density (\(N_{\rm S}<1\)) will show the images of individual scatterers, while images with a high source density (\(N_{\rm S}>1\)) will show speckle, as the signal from the individual scatterers will interfere. Due to the relatively low resolution in echography speckle will likely be prominent in many UIV studies, unless very low tracer/scatterer contractions are used. Increasing the frequency will reduce the wavelength (and thus improve the resolution), but also strongly increase the source density.

Despite the fact that speckle is the result of interference of many particles, it is still possible to estimate the local displacement from the ensemble of scatterers, as long as the displacements are small. Larger displacements and gradients in the flow will make the matching of patterns more difficult: a small change in the configuration of the scatterers can lead to a significantly different speckle pattern (Friemel et al. 1998). For images in the low source density regime (\(N_{\rm s}<1\)), this is generally much less of an issue. Nevertheless, in both regimes the same principles can be used for pattern tracking to find the local displacement and thus velocity.

As mentioned, speckle tracking has been used to quantify the motion and deformation of tissue (e.g. muscle contraction); this approach is generally referred to as ‘block matching’. A small subregion of the B-mode image is selected and an algorithm then finds the region in the subsequent frame that has the best match with the original gray-scale pattern. Various approaches to achieve this exist, including minimising the summed error, optical flow, and correlation functions Wang and Shung (1996), Hein and O’brien (1993). The latter approach is the most common and illustrated in Fig. 12. A region, known as either an interrogation area, window, kernel or block, is extracted from two frames. The cross-correlation of these two areas is referred to as the correlation plane. For a purely uniform displacement this correlation plane will contain a dominant peak, located at dx and dy from the origin, with the coordinates representing the most likely displacement between the two. When an error-minimisation strategy would have been used, this displacement would result in a minimum in the penalty function (e.g. based on the summed difference in intensity). Dividing these displacements by the time between the frames (\(\Delta T\)) and taking into account the calibration/scaling factor then gives the local velocity vector. This process is repeated for all other regions in the original frame, so that an instantaneous vector field is obtained. By utilising the finite width of the correlation peak, it is possible to estimate the displacement with sub-pixel accuracy by using some sort of peak fit. The interrogation area requires a certain size, so that it contains sufficient image features that can be tracked. The optimal size depends on, among others, the scatterer concentration, the nature of the flow, and the image quality. In practice windows of 8 × 8 to 64 × 64 are generally used. This means that the resolution of the velocity field is generally an order of magnitude lower than than the resolution of the images. Note that this is in contrast with e.g. MRI, where each voxel provides a velocity estimate. For an in-depth discussion of PIV, including more sophisticated algorithms, the reader is referred to e.g. Adrian and Westerweel (2010).

Fig. 12
figure 12

Reproduced from Poelma and Fraser (2013). ©IOP Publishing. Reproduced with permission. All rights reserved

Schematic representation of particle image velocimetry processing of a B-mode image pair (shown left as two superimposed frames): the local cross-correlation result (middle) provides the most likely displacement within each interrogation area (IA). The final result is an instantaneous vector field (right)

Various improvements for speckle tracking were introduced over years, such as the use of multi-level block sizes (Yeung et al. 1998) and deforming meshes (Yeung et al. 1998); these are equivalent to the development of adaptive interrogation area sizes and window deformation in PIV (Adrian and Westerweel 2010). As stated above, the correlation approach works for both speckle and low source density images, where individual scatterers can be observed. The main difference will be the more rapid ‘decorrelation’ for the speckle images, so that the displacement and thus \(\Delta T\) requirements are more stringent (Adrian and Westerweel 2010; Friemel et al. 1998). A consequence is that adding more tracer (e.g. microbubbles) may actually decrease the signal-to-noise ratio beyond a certain threshold (Niu et al. 2011), as the nature of the image changes to a high source density image. Note that some studies actually make use of the rate of decorrelation to estimate the velocity magnitude (Bamber et al. 1988; Rubin et al. 1999).

The first studies reporting blood flow estimation using cross-correlation of B-mode images were by Trahey et al. (1987, 1988) at Duke University in the late 1980s, summarised in the review by Bohs et al. (2000). Correlation-based techniques were introduced specifically to create an ‘angle-independent’ method, i.e. a method that yields two velocity components, regardless of the angle between the acoustic beam and the flow. A result from Trahey’s experiments is reproduced in Fig. 13; this particular figure only shows a profile, but the study also reported vector fields. Note that this first ‘in vivo PIV using ultrasound’ study was only a few years after the first application of PIV in fluid mechanics. The correlation-based approach by Trahey and coworkers did not find widespread use, as it was computationally intensive Wang and Shung (1996) and thus not deemed suitable for real-time clinical applications.

Fig. 13
figure 13

Reproduced from Trahey et al. (1987). ©1987 IEEE. Reprinted, with permission, from IEEE Transactions on Biomedical Engineering

An early example of blood flow measurement using correlation of B-mode images. The image shows a human popliteal vein, acquired using a 10 MHz scanner. The transducer position is located at the top of the image. The vectors represent the result of the local correlation search

The main breakthrough of UIVFootnote 6 came in the early 2000s. By that time, PIV had become commonplace in fluid mechanics laboratories. A sound understanding of the underlying theory and many new algorithms had matured it in an accurate and flexible tool. Researchers started applying standard PIV algorithms on data obtained from echography, realising that this would enable them to measure in opaque fluids. One of the first reported examples was a study in sediment transport by Crapper et al. (2000), who used a standard medical system at 3.5 MHz and VidPIV software. Clear water experiments were performed using conifer pollen as tracer, while in the mud suspensions the dispersed phase (16 g/l) functioned as tracer. Their ultrasound-based results for the clear flows matched within 15% of reference measurements using conventional PIV. For higher velocities, they noted that the tracer images became more streak-like.

The development of an ultrafast ultrasound system, originally intended for elastography, inspired Sandrin et al. to measure speckle displacement in a corn starch solution. Using a frame-rate of up to 5000 1/s, they obtained the displacement of speckle in the axial (beam) direction by means of cross-correlation. This made it possible to obtain snapshots of the velocity field of one component in a stretched vortex (Sandrin et al. 2001). In a follow-up paper, they applied two-dimensional cross-correlation on the RF data to obtain both velocity components (Manneville et al. 2001). They also obtained results using whole milk as working fluid (Sandrin et al. 2001; Manneville et al. 2001) and using a suspension containing Levovist contrast medium (Manneville et al. 2001).

Around the same time, Gharib and coworkers at the California Institute of Technology started applying ‘digital ultrasound speckle imaging velocimetry (DUSIV)’ to flow phantoms representing the aorta (Zarandi et al. 2001; Lin et al. 2003; Gharib and Beizaie 2003). The use of echography (based on a commercial medical ultrasound system) allowed them to quantify the flow through non-transparent, compliant walls of the model. No tracer was needed, as the system used speckle from the bovine blood that was used as medium.

The biggest contributor to the popularity of UIV was the seminal paper by Kim et al. (2004a) at the University of Colorado. This paper clearly laid out the basics and showed the potential of ‘echo-PIV’ for the fluid mechanics community. It demonstrated that it could be accurate, as a comparison with optical PIV measurements and theoretical profiles in a fully-developed laminar flow showed an agreement within 7.7%. Images were obtained again using a medical ultrasound system operating at 10 MHz and the study emphasised the value of tracers specifically for UIV (microbubbles, see next section). A number of follow-up papers (e.g. Kim et al. 2004b; Liu et al. 2008) from the same group explored more advanced implementations, such as the possibility to obtain real-time measurements (Zheng et al. 2006). An example taken from the latter study is shown in Fig. 14, where a vortical flow is measured using both optical and echo-PIV.

Fig. 14
figure 14

Image reproduced from Zheng et al. (2006), with the permission of AIP Publishing. ©2006 American Institute of Physics

Validation of echo-PIV using a vortical flow: (a) B-mode particle image of the flow; (b) velocity field measured by echo-PIV; and (c) echo-PIV and optical PIV velocities along one radial line within the flow field

The paper by Kim et al. inspired several other groups to apply UIV in a range of applications, a number of these are summarised in Sect. 5. Most of these studies used conventional PIV software on data obtained using standard medical ultrasound devices. While this has been a successful strategy, there are some distinct differences and potential pitfalls that must be avoided. These are discussed in the subsequent section.

4 Major difference between UIV and PIV

4.1 Sweep correction

A fundamental difference between PIV and UIV is the timing of the image formation. In PIV, cameras obtain a snapshot of the tracer images. In echography—with the exception of plane-wave imaging—the image is constructed by sequentially recording scan lines. This implies that different parts of the image are recorded at different times, similar to a rolling shutter. While in PIV the velocity is obtained simply by dividing the displacement by the delay time between frames, in UIV the delay time is actually also a function of the displacement. This becomes more and more apparent for objects moving with a lateral velocity that approaches the sweep velocity, the hypothetical ‘velocity’ at which the system progresses to the next scan line. This sweep velocity can be determined by multiplying the lateral size of the domain with the frame rate.

The consequences of the ‘beam sweep’ effect were acknowledged by some authors (Beulen et al. 2010; Poelma et al. 2011), and an extensive discussion is given in the study by Zhou et al. (2013). The velocities can be corrected for the beam sweep effect using the following expressions:

$$\begin{aligned} V_x = \frac{V_{\text{s}} V_x^{\prime}}{V_{\text{s}}+V_x^{\prime}},\ \ \ \ V_y = \frac{V_{\text{s}} V_y^{\prime}}{V_{\text{s}} + V_x^{\prime}} \end{aligned}$$
(5)

Here \(V_x\) and \(V_y\) are the corrected velocity components, \(V_x\) and \(V_y\) the uncorrected values, and \(V_{\text{s}}\) is the sweep velocity. This correction is exact for flows without acceleration. Note that an axial velocity (\(V_y\)) does not lead to a time difference by itself (it is absent in the denominator), but a correction is still needed to accurately estimate this velocity component. Furthermore, it is important to take into account the sign of the velocity, i.e. whether the tracers are moving with or against the sweep direction. In most cases the recorded image data contains no information about the sweep direction. It should thus be confirmed at the time of acquisition. Some transducers have an orientation marker, but a simple alternative is to record the same flow using a flipped transducer and evaluate the difference in mean displacement.

In practice, the relative errors in the velocity can be up to 20% without correction for low frame rates and high velocities (Zhou et al. 2013). Zhou et al. also predicted that for transient flows additional corrections are needed. Recently, Fraser et al. (2016) showed that in pulsatile flow in a rabbit aorta this correction was indeed essential to obtain agreement with reference values. It is unclear whether the studies that do not mention the sweep effect perform a correction for it or not. It may explain some discrepancies in validation studies if it was omitted. The sweep effect may also explain the observed ‘streaky’ particle image shapes reported by Crapper for high velocities (Crapper et al. 2000). The correction method reported in Eq. 5 only holds for conventional image formation. The processing of images obtained using e.g. parallel beam forming (where two or more scan lines are recorded simultaneously across different parts of the field-of-view) and interleaved imaging requires specific corrections. Derivation of these are trivial, as long as the exact time stamp for each scan line is known.

4.2 Non-uniform image resolution

In conventional PIV images are created by aligning the focal plane of the optical system with the plane formed by the thin laser light sheet. As only particles within the thin light sheet are illuminated, all tracer particles appear focussed in the image captured by the camera. Under normal conditions, the spatial resolution is constant throughout the image. For echography, this is no longer the case. Generally, a single focus is placed at the centre of the field-of-view. As shown schematically in Fig. 8, tracer particles that are located before or beyond this focal location (where the resolution is optimal, i.e. the Azimuth Beam Width, ‘ABW’) will have a larger image than those in the focal region. An example of this effect is shown in Fig. 15, where the point-spread-function (PSF, the observed image of a point source) is shown as a function of distance. Near the focal plane the resolution is optimal, with the smallest PSF. Closer to the transducer the PSF has an elongated shape: the axial resolution is constant throughout the field-of-view, while it become worse in the lateral direction away from the focal point. Around an axial position of 20 mm the lateral resolution is approximately 4 times lower than at the focal location (28 mm). The non-uniform resolution is also evident in e.g. Fig. 4, where the region near the transducer appears somewhat more blurry than the region in the middle of the image. This effect is distinct from the inherent non-uniform resolution that occurs when scan conversion is applied to correct for the transducer shape.

Fig. 15
figure 15

Image reproduced from Contreras Ortiz et al. (2012). ©2012 Elsevier Ltd. Used with permission

Variability of the Point Spread Function in the axial direction. A phantom is constructed with 0.1-mm-diameter nylon wires. This phantom is imaged using a 7.5 MHz linear transducer, with its focus placed at an axial position of 27 mm

While it may be tempting to use more focal regions, this should be avoided for velocimetry. First of all, the imaging rate reduces linearly with the number of focal regions, as different parts of the image are acquired sequentially, before they are merged. Apart from the frame rate issue, the sequential process also introduces a time lag between different regions of the image. This effect is shown in movie S1 (supplemental material), where particles can be seen ‘hopping’ from one part to the other (at the positions noted by the triangles), as a result of the time delay between acquisition. This will introduce artefacts in the velocity field.

Apart from the non-uniform in-plane resolution, the variation in the elevational width of the field-of-view should also be considered. While normal PIV has a well-defined, thin measurement volume, in echography the focussing in the elevational direction is done using a fixed acoustic lens (see also Fig. 2). At the focal position the width is typically a few millimeters, but it steadily increases away from the transducer. This means that velocity data is obtained as an average over a fairly thick slice. Very few papers report the actual thickness of the slice, while it is a critical parameter when interpreting results.

A final consideration is the inherent difference in resolution in the lateral and beam directions. The resolution in the lateral direction is determined by the transducer pitch and focussing details (as discussed earlier). The resolution in the beam direction is determined by the wavelength and pulse details. It is generally also sampled at a higher frequency than the original ultrasound frequency. This means that the apparent resolution can be an order of magnitude different; this implies that the aspect ratio of a pixel is far from unity. If this data is processed in this form by PIV software, tracer particle images will appear stretched. This may lead to issues with the sub-pixel displacement estimators of conventional PIV algorithms: these are generally built on the assumption that particle images are ‘small’ and Gaussian in shape. Depending on the entire imaging chain (acquisition, line density, beamforming, interpolation) the actual data may be far from this ideal situation, as is also evident from Fig. 15. This can introduce severe pixel locking effects Adrian and Westerweel (2010). This effect, which results in a strong bias towards integer displacements, leads to distinct ‘steps’ in velocity data, see also Fig. 14c. It is possible to reduce these effects to a great extent by using alternative sub-pixel displacement estimators, optimised for the nature of the data (Poelma et al. 2012). Another consequence of the pixel aspect ratio is that the same physical displacement in the beam and lateral direction will lead to a pixel displacement that can be different by an order of magnitude (Gurung and Poelma 2016). Again, this may lead to problems with conventional software (e.g. during iterative processing with window deformation or during the vector validation stage). Both these problems can be alleviated by resampling the data at the actual resolutions, whose ratio is generally much closer to unity. As an additional bonus, this will reduce the processing time needed.

Only a limited number of studies have investigated the differences between using RF data or data that has been processed to ‘B-mode’ images (Chen et al. 2010). Generally, steps such as envelop detection and log-compression of the image have little effect on the displacement estimation Gurung and Poelma 2016). This can be understood by realising that these steps may alter the shape of features in an image, yet not their location and do so in both frames identically. As long as the processing steps retain traceable image features, only minor differences are to be expected, likely due to the performance of the sub-pixel estimation. For most studies these effects will not be evident in e.g. mean flow statistics. Nevertheless, more extensive investigation of the pre-conditioning of data for velocimetry may be fruitful in the future (Yeom et al. 2014a). Naturally, image enhancement options like Speckle Reduction and Clutter Rejection or lossy image compression should be avoided at all cost.

4.3 Seeding material

Tracers for UIV need to satisfy the same basic requirements as conventional seeding material for PIV: they need to scatter efficiently and should faithfully follow the flow. The latter requires that the Stokes response time of the particleFootnote 7 is sufficiently small to follow even the fastest fluctuations in the flow and there should be no appreciable settling (or rise) velocity. As all UIV experiments have a liquid as continuous phase and are generally in a low Reynolds number regime, the response time criterium is easily met. Naturally, the tracers should not affect the fluid phase—unless this is the purpose of the study, as discussed later.

As discussed in Sect. 2.3, the scattering efficiency is proportional to the change in acoustic impedance. Several particle types commonly used in PIV have been used successfully in UIV, such as Sphericel (Gallot et al. 2013; Gibaud et al. 2016) and Vestosint Poelma and Gurung (2016). Much better signal-to-noise ratios can be obtained by making use of tracer materials specifically developed for ultrasound imaging: ultrasound contrast agents (UCAs). These are micrometer-sized bubbles of air or heavy gases such as octafluoropropane (\(\hbox {C}_3\hbox {F}_8\)) or sulphur hexafluoride (\(\hbox {SF}_6\)). A lipid shell enhance the lifetime of the bubble, which would otherwise dissolve in seconds (Jong et al. 2009). Their acoustic impedance leads to strong scattering, which is enhanced by non-linear oscillations of the bubble shape (Frinking et al. 2000; Mukdadi et al. 2004). UCAs are commercially available through various companies and various types/brands have been used in UIV experiments. Examples include Optison (GE Healthcare, Kim et al. 2004a; Zheng et al. 2006; Westerdale et al. 2011), Definity (Lantheus Medical Imaging, North Billerica, MA; Kheradvar et al. 2010; Walker et al. 2014), Levovist (Schering/Bayer; Manneville et al. 2001) and SonoVue (Bracco; Zhang et al. 2011; Poelma and Fraser 2013; Poelma et al. 2011). As these UCAs are approved for use in patients, their manufacturing process is very strictly controlled and thus they are fairly expensive (up to hundreds of euros for a 2 ml vial). However, it is also possible to make them in-house if larger quantities are needed (Qian et al. 2010). Note that UCAs in UIV applications are used in a much lower concentration than in its original intended use: for contrast-enhanced ultrasound imaging, typically a highly-concentrated bolus is injected to get an optimal visualisation of blood flow. Using a more diluted, continuous feed (infusion) it is possible to obtain much better data quality for UIV purposes (Cimino et al. 2012; Fraser et al. 2016).

The effect of the concentration of contrast microbubbles was studied by Kim et al. (2004a). However, as the concentration of the bubbles was strongly time-dependent, it was difficult to formulate straightforward guidelines; they found the lowest error at concentrations around 103 bubbles per mL. Based on the data available for that study, this corresponds to an estimated source density \(N_{\rm S}\) of order unity, denoting the transition to speckle (see Sect. 3.3). Using synthetic data Gao et al. (2011) also investigated the effect of the concentration. For the simulated flow in a ventricle, they obtained the best results with 10 bubbles per mL. Niu et al. found an optimal concentration for an in vitro experiment of 0.8–2 × 103 and 1–5 × 105 bubbles/mL for an in vivo experiment. As the optimal concentration is clearly different for various implementations, it is advisable to establish it for each new case. A real-time method to achieve this was introduced by Niu et al. (2011) based on a number of image characteristics. Recently, Walker et al. found that the concentration of contrast microbubbles was not uniform in a flow phantom of a stenosis. They speculated that the flow-bubble interaction excluded them from a recirculation zone, despite the low Stokes number. As the exact physical details of the bubbles are generally not available, it is difficult to predict the response of microbubbles to strong flow gradients.

For applications with a heterogeneous fluid—e.g. multiphase systems that inherently contain scatterers—no seeding needs to be added. For instance it is possible to directly track speckle of red blood cells (Nam et al. 2012; Yeom et al. 2014b; Azevedo et al. 2015), see also Sect. 5 and the work by Trahey et al. discussed earlier. Similarly, suspensions such as milk (Sandrin et al. 2001) and dense concentrations of e.g. clay have sufficient inherent scattering behaviour to obtain velocity fields (Gurung et al. 2016), see e.g. Fig. 1. For applications such as sediment transfer, the dispersed phase is generally used to determine the velocity field (Crapper et al. 2000; Chinaud et al. 2010). As the Stokes number will be larger than unity, the measurements will no longer represent the (carrier) fluid velocity, but only the dispersed phase. This can make the comparison of single-phase and particle-laden cases difficult (Poelma and Gurung 2016).

To conclude the discussion of tracer material, it is worthwhile to briefly mention some common imaging artefacts that may occur. An extensive discussions can be found in the reviews of e.g. Contreras Ortiz et al. (2012) or Feldman et al. (2009). For UIV, it is desirable to obtain small, distinct tracer images, such as shown in Fig. 14a. When a relatively large object (i.e. larger than the resolution) consisting of an attenuating material is present in the field-of-view, it can create a shadow. An example is shown in Fig. 16a, where small air bubbles are present in the image. The opposite can also occur, when the object is of a material that is less attenuating, it creates a lighter band [increased through transmission Feldman et al. (2009)]. In Fig. 16b reverberation effects can be seen. In this case, the flow contains poorly-wetted polystyrene particles; small air bubbles remain attached to their surfaces. Rather than creating a single reflection, multiple echoes are created, possibly reinforced by bubble resonance (ring down artefact). These reverberations occur after the original reflection, so that they appear below the particles, creating a comet tail. The effective axial resolution will become worse as a result of these comet tails. When the particles are completely wetted by means of a surfactant these artefacts disappear, as is evident in Fig. 16c. Similar reverberations, caused by sound reflecting multiple times are also often evident near walls, see e.g. Fig. 3d.

Fig. 16
figure 16

Imaging artefacts: (a) shadowing due to medium-sized air bubbles; (b) comet tail/reverberations from poorly-wetted polystyrene particles; (c) the same polystyrene particles after proper wetting

5 Application examples

In this section, a selection of applications of UIV is presented to show the versatility of the technique. Section 5.1 discusses a number of validation studies in various flow geometries. Examples of cardiovascular flow applications are given in Sects. 5.2 and 5.3 discusses some more conventional engineering applications.

5.1 Validation studies

Understandably, the first wave of studies using UIV focussed on the validation of the technique. This was done for increasingly complex flows: initially using steady laminar flow (Beulen et al. 2010; Zheng et al. 2006; Qian et al. 2010; Niu et al. 2010; Poelma et al. 2012), followed by pulsatile flows (Kim et al. 2004b; Poelma et al. 2011) and ultimately in flow phantoms. The latter are ‘in vitro’ studies using models that represent e.g. a ventricle (Kheradvar et al. 2010) or (diseased) blood vessels (Zhang et al. 2011; Zhu et al. 2011). Validation was done by comparing the UIV velocity fields with theoretical profiles (for the case of Poiseuille flow), by comparing with CFD results (Liu et al. 2008; Qian et al. 2010; Zhu et al. 2011) or by direct comparison with conventional PIV (Kim et al. 2004b; Zhang et al. 2011; Westerdale et al. 2011; Kheradvar et al. 2010; Walker et al. 2014).

Generally, the agreement between the results obtained with UIV and reference data was found to be satisfactory, with discrepancies well below 10%. On average an error of approximately 5% was reported. Note that these low discrepancies only represent optimal situations, where the acquisition rate was sufficiently high to capture the flow. Under well-controlled conditions and using careful analysis (i.e. including sweep correction), random errors in the velocity can be reduced to 1% (Poelma et al. 2012). This shows that, in principle, UIV can achieve an accuracy comparable to conventional flow measurement techniques. However, the resolution cannot match that of e.g. conventional PIV. Note that most of these studies rely on (phase-locked) averaging. This averaging can be done either based on the velocity vectors or the intermediate correlation planes (Poelma et al. 2011; Chen et al. 2010). Nevertheless, recent work has shown that also instantaneous measurements can be accurate within a few percent: Gurung and Poelma (2016) compared the statics of a fully-developed turbulent pipe flow with literature data and found a good agreement.

5.2 Cardiovascular flows

The vast majority of UIV applications studies reported in the literature focusses on cardiovascular flows, which can easily be explained by the fact that the required hardware and expertise is readily available in this field. Studies mostly focus on two main goals: (1) characterise the hemodynamic forces in blood vessels to study the role of blood flow in the inception and progression of diseases (2) study the hemodynamics inside the heart.

Blood flow inside the heart has been studied by a collaboration between Pedrizzetti, Sengupta and Kheradvar (and other co-authors) in a series of papers (Sengupta et al. 2007; Hong et al. 2008; Kheradvar et al. 2010; Cimino et al. 2012). After validation using an in vitro model, they reported the ventricle filling process in healthy volunteers and diseased patients. Using multiple planes (either parallel or orthogonal) they were able to reconstruct the three-dimensional flow field (Sengupta et al. 2012; Pedrizzetti et al. 2014). To better understand the interaction between blood flow patterns and the moving wall, they also determine the pressure field from the velocity data (Cimino et al. 2012). UIV is uniquely capable of providing out-of-plane vorticity data, as two components of the velocity field are recorded. This makes it an excellent tool to study vortex dynamics in the ventricle, which is thought to be an indicator of cardiac health (Pedrizzetti et al. 2014). A typical result is shown in Fig. 17. This figure shows the flow pattern in the left ventricle of the heart of a healthy volunteer during a cardiac cycle (from ejection to filling). The data was obtained using injection of UCAs and imaging in the 3–5 MHz range at 80 Hz.

Fig. 17
figure 17

Image reprinted from Hong et al. (2008) with permission from Elsevier. ©2008 American College of Cardiology Foundation

Time sequence of the ejection (ac) and filling during the isovolumic relaxation of the left ventricle (e, f). The vectors represent blood flow velocities and are superimposed on B-mode imageas. The wall is indicated by the purple line

Studies of blood flow have been reported in rats (Qian et al. 2010) and in rabbits (Qian et al. 2010; Fraser et al. 2016). So far these were primarily proof-of-principle experiments to show that the mean flow statistics match with reference (Doppler) data. The motivation of using UIV is the fact that it can provide the wall shear stress in complex geometries (Poelma et al. 2012), which is an important parameter in e.g. the onset and progression of cardiovascular diseases such as atherosclerosis and aneurysms.

A final example of an application of UIV in cardiovascular flows is shown in Fig. 18. This figure, taken from the study by Nam et al. (2012) shows the pulsatile flow through the valves in a vein near the ankle of a human volunteer. These in vivo snapshots show the wealth of data that is available for biomechanical fluid-structure interaction studies: apart from the flow field, the position of the valves can clearly be observed in the B-mode images. The data was obtained using high-frequency ultrasound (35 MHz), so that speckle of red blood cells could be used for PIV analysis. Note that the vein has a diameter of approximately 2.5 mm, indicating that this experiment had a much smaller field-of-view than previous examples. It should be pointed out that ‘fast’ parts of the cardiac cycle could not be obtained due to the low imaging rate (which would be inbetween panels c and d). With an interframe time of 33 ms, only velocities up to 4–5 mm/s could be obtained.

Fig. 18
figure 18

Reproduced from Nam et al. (2012). ©Springer Science+Business Media, B.V. 2010. Used with permission

Instantaneous velocity fields of the perivalvular blood flow superimposed on B-mode images of 500 × 500 pixels during the opening and closing phases of the valve in the greater saphenous vein at the ankle; ac correspond to the valve opening stage, while df depict the valve closing stage

5.3 Miscellaneous applications

The earliest example of UIV (apart from Trahey’s work) was in the field of sediment transport (Crapper et al. 2000). This application is an evident candidate for UIV: dense concentrations of particles make optical techniques impossible, while the particles can act as tracer for ultrasound-based methods. Crapper et al. was able to obtain velocity fields in densely-laden suspensions flowing past a baffle in a flume. It should be stressed once more that this approach provides the velocity field of the dispersed phase (solid particles) and not the fluid.

Only recently other researchers followed the approach of applying ultrasound imaging in this application area. Zou et al. (2015) characterised the topography of a sediment bed and also used the technique to study the onset of sand motion (‘sediment inception velocity’) (Song and Zou 2015). In another study, they used the local scattering intensity (i.e. the brightness in the B-mode image) to determine the local sediment concentration (Zou et al. 2014). Calibrations where performed with uniform suspensions of known concentration, which provided a calibration curve for concentration versus brightness. A similar approach was followed by Gurung and Poelma (2016), who found that the intensity (I) could be described by a function of the form \(I \sim 1-{\text{exp}}(-c/c_0)\), with the intensity reaching a saturation value around a volume fraction of 10%. They obtained concentration profiles and mean velocity profiles in horizontal, particle-laden pipe flows for two cases: the case were a mobile bed was present and the case where all particles remain suspended. Figure 19 show two B-mode images of both cases, together with the measured mean velocity profiles. Turbulence statistics could be obtained in this flow at moderate volume fractions up to at least 3%.

Fig. 19
figure 19

Reproduced from Gurung and Poelma (2016) under the Creative Commons Attribution 4.0 International License. ©A. Gurung and C. Poelma

Exemplary raw images obtained in a flow with a partially mobile particle layer (left) and with fully suspended particles. For both cases the mean streamwise velocity profile is also shown

Another area where UIV can be useful is in the study of non-transparent, complex fluids, such as foams, suspensions and colloidal systems. Using a combination of a velocity profile measurement and a pressure drop measurement, it is possible to determine the rheological behavior of these fluids. Early implementations of this method used Doppler profiling [‘ultrasound Doppler velocimetry pressure drop’, UDV-PD method (Wiklund et al. 2007; Ouriev and Windhab 2002)]. Once a velocity profile is obtained and a constitutive equation for the fluid is assumed, the rheological parameters follow from a fit of the analytical solution for the velocity profile and the known pressure drop. This approach can be used in cylindrical pipes for in-line rheology applications. Examples include e.g. studies of paper pulp up to a volume fraction of 7.7% (Wiklund et al. 2006), industrial suspensions and foods (Wiklund and Stading 2008; Wiklund et al. 2010). These types of flow are echogenic by nature, i.e. they contain sufficient acoustic scatterers so that no seeding needs to be added. Recently, Gurung et al. (2016) showed that UIV (rather than Doppler) can be used in similar rheology studies. An example is shown in Fig. 1, where the velocity field is measured in a model drilling fluid flowing steadily through a pipe. While the flow in this figure only has a streamwise (lateral in the measurement frame) component, UIV can also characterise, for instance, the flow encountering a sudden diameter change (Gurung et al. 2015). This opens up more sophisticated methods to establish rheological parameters (Drost and Westerweel 2013).

Ultrasound-based velocimetry has also been successfully applied in Taylor-Couette flows, for more fundamental rheology studies (Manneville et al. 2004). In this case the torque, rather than the pressure drop, is used to quantify viscous losses. Ultrasound speckle velocimetry was used to characterise fatigue due to shear in gels (Gibaud et al. 2016), flow banding in cellulose nanofibril suspensions (Martoïa et al. 2015), flow instabilities and shear banding in a non-Newtonian micellar solution (Gallot et al. 2013) and avalanche-like fluidisation events in suspensions (Kurokawa et al. 2015). Again, most of these systems are inherently echogenic and do thus not require seeding. A drawback of the Doppler (and 1-D speckle) approach, in particular for the Taylor-Couette geometry, remains the fact that only the component in the beam direction can be obtained. Apart from making the technique less sensitive (as a projection of the streamwise or tangential velocity is measured), it prohibits the use of this approach in e.g. turbulent flows. Nevertheless, in all of these studies the velocity fields provided insight in the underlying physics of the observed overall rheological behaviour, hereby “linking global rheology to the local dynamics” (Kurokawa et al. 2015). This makes ultrasound-based techniques a very powerful tool in rheology.

6 Conclusion and outlook

Ultrasound imaging velocimetry (or echo-PIV) has, over the last decade, matured to a state where it can provide valuable insight in flows that are not accessible for optical measurement techniques. It can provide accurate, instantaneous two-component velocity fields with a resolution of a few millimeters, across a field-of-view on the order of several centimeters squared.

Initially, measurements were performed by directly applying conventional PIV software on image data exported from existing medical ultrasound devices. Recent work has shown that better results can be expected when all steps of the process are optimised for UIV. This includes: (1) designing dedicated test sections for acoustic access, rather than optical access; (2) using specific tracers, such as contrast agents; (3) usage of unprocessed (RF) data; (4) image acquisition hardware and/or protocols to enable short inter-frame times; (5) cross-correlation algorithms that are tuned for echography data, including e.g. sweep correction and robust sub-pixel estimators. All of these steps are topics of current research and significant improvements can be expected in each step individually, but also in their integration.

A major opportunity exists for better utilising the raw RF data. This data is generally converted to an image and processed, while other processing methods can be formulated specifically for data of this nature (Bodnariuc et al. 2015). The original RF signal contains more information than the B-mode images. This information can be used to precondition the data (Yeom et al. 2014a), to address some of the inherent ultrasound imaging issues highlighted in Sect. 4.2. Alternatively, this information can be used to extract specific details. In the original paper by Kim et al. (2004a), the possibility of using harmonic imaging was already explored. Contrast microbubbles non-linearly oscillate in response to ultrasound of certain frequencies. By isolating the second harmonic in the recorded data, Kim et al. could enhance the signal from the microbubbles and suppress other signals (e.g. reflections at the original frequency). Apart from a better signal-to-noise ratio, this strategy has also been suggested to measure both the fluid and dispersed phase in particle-laden flows (Gurung and Poelma 2016), similar to the use of fluorescent particles in optical PIV.

Two recent major technological developments in echography will benefit UIV: plane-wave imaging and matrix transducers. Plane-wave imaging allows the acquisition of images with frame rates in the kHz range. UIV benefits from this development, as it enables measurement in faster flows (Rodriguez et al. 2013; Gallot et al. 2013; Leow et al. 2015). Furthermore, it enables the use of sliding correlation averaging in temporally oversampled data. As the signal-to-noise ratio in UIV is generally rather low, averaging results over a number of frames will drastically improve the data quality in transient flows (Poelma et al. 2011). Matrix transducers, i.e. transducers with a grid of elements, have been introduced to obtain 3D/4D image data. In combination with UIV this allows 3D flow measurement without the need to scan plane by plane. A drawback is the even lower frame rate compared to linear transducers. Furthermore, the number of elements is usually limited to avoid excessively expensive electronics. A feasibility study based on synthetic data was reported by Gao et al. (2013) and the first demonstration in practice has also been presented (Falahatpisheh and Kheradvar 2014).

Currently, most applications focus on cardiovascular flows. This is understandable, considering the history of the technique. However, there are a myriad of other applications that can benefit from UIV. These range from studies of sediment transport to flows of cosmetics, food products, or other complex fluids. The studies summarised in this review show that UIV can be an accessible, reliable flow measurement technique; it is therefore likely only a question of time before the technique finds its way to a wider audience.

7 Supplemental data

7.1 Movie S1

This movie shows the acoustic streaming effect: by using a high powered ultrasound signal (at 10 MHz), a flow is generated, pushing the fluid and Vestosint tracers away from the transducer (located outside the frame, on the top). The field-of-view is approximately 4 × 4 cm2 and the maximum velocity is approximately 1 cm/s. The data was acquired at 31 fps.

Also shown is the effect of using three focal points while acquiring images. The transition from the merged regions is indicated by the yellow triangle. Particles can be seen to ‘hop’ from one region to the other, as the three focal regions are not recorded simultaneously, but sequentially.