A Thermoviscous Analysis of Acoustic Radiation Forces

This past July, I had the pleasure of attending the 22^nd International Congress on Sound and Vibration. In addition to running the COMSOL vendor booth with my Italian colleague Gabriele, I was also a presenter at the event. My presentation was based on a paper I wrote with Henrik Bruus and Jonas Karlsen that focuses on how to determine acoustic radiation forces including thermoviscous effects. Let’s explore acoustophoretic effects in greater detail and the research findings highlighted in my presentation.

What Is Acoustic Radiation Force?

Acoustic radiation force arises due to momentum transfer from an acoustic field to a small particle, which results in a net force acting on the particle. In the related phenomenon of acoustic streaming, a net fluid flow is driven by an acoustic field, with momentum transfer (to the fluid) occurring in the acoustic boundary layer.

Both phenomena are important as they are used to move, control, and sort particles or small living organisms (i.e., cells, proteins, or large viruses) in, for instance, microfluidic devices. This kind of particle manipulation is known as acoustophoresis, that is, movement by sound. The radiation force acts directly on the particle while acoustic streaming induces drag on the particles. The magnitude and/or direction of the effects depends on the size and mechanical properties of the particles as well as the fluid properties. A solid understanding of such phenomena is essential for producing and optimizing devices to efficiently sort and manipulate particles. In a previous blog post, I demonstrated how to model the acoustic radiation force and acoustic streaming in a microfluidic device.

My presentation at the 22^nd International Congress on Sound and Vibration (ICSV22) was an extension and generalization of how to model the acoustic radiation force in detail using the so-called second-order perturbation theory. The Acoustophorectic Force blog series provides a helpful introduction to understanding the radiation forces as well as different ways to compute these forces.

Different Approaches to Analyzing Acoustic Radiation Forces

Let’s now get a bit more technical and pull out a few equations. As is the case for most physics problems, there are several ways of describing and deriving the radiation force effect under different physical assumptions and in different regimes. I will briefly summarize and describe these different approaches here.

Direct FSI Approach

In the direct fluid-structure interaction (FSI) approach, the full (nonlinear) set of governing equations is solved as a dynamic problem. The radiation force is then extracted by looking at the time-averaged particle motion. In order to include both thermal and viscous losses, the Navier-Stokes equations need to be solved together with the energy equation, which is, in essence, nonisothermal flow. The equations should be solved with high numerical precision in order to capture all of the small-scale effects associated with the radiation force.

Second-Order Perturbation Approach

To simplify the analysis of the full set of nonlinear governing equations, the problem can be split up into a linear acoustic problem and a linear time-averaged second-order problem. This is accomplished via the second-order perturbation theory; for example, the pressure p is written as p = p_0+p_1+p_2. Subscript 0 represents background properties (the 0^th pressure could be roughly a constant single atmosphere, depending on equilibrium conditions), subscript 1 represents the classical acoustic fields, and subscript 2 represents the second-order perturbation fields needed to capture the nonlinear effects. The different-order equations are solved sequentially.

This theory establishes that, in general, the radiation force on an object \Omega is given by

(1)

where \langle \cdot \rangle represents time averaging, \sigma_2 is the second-order contribution to the stress, \mathbf{u}_1 is the acoustic variation in the velocity field, and the integral is on any closed surface surrounding the particle.

Thus, in order to find the radiation force, the acoustic problem (first-order perturbation) and the second-order perturbation equations must both be solved. This method was used in the research presented at ICSV22, which we will discuss in greater detail below. The equations are solved, including thermal and viscous losses in both fluids and solids.

The No Loss Limit

The fluid can be treated as having no loss for particles of size a much larger than the viscous and thermal boundary layers (the penetration depth \delta_\textrm{tv}), that is, \delta_\textrm{tv}/a \ll 1. In this case, the second order-stress in Eq. (1) can, in the far-field, be simplified to

(2)

By determining the scattered field p_1 and \mathbf{u}_1 of any object, the radiation force in Eq. (1) can be calculated. It is not necessary to know the second-order fields explicitly. Ref. 2 in our Further Reading section and the blog post “How to Compute the Acoustic Radiation Force” provide two examples of this.

The Rayleigh Limit

In the Rayleigh limit, the particle is assumed to be acoustically small. That is, k\cdot a \ll 1, where a is the particle size and k is the wave number. Scattering theory is used to derive the scattered field off the particle. In this limit, the radiation force given in Eq. (1) reduces to

(3)

where f_0 and f_1 are the monopole and dipole scattering coefficients and the subscript \textrm{in} represents the incident acoustic field (e.g., the standing acoustic field in a sorting device) evaluated at the small particle’s location. The fluid compressibility is \kappa_s and the fluid density is \rho_0.

The whole gymnastics of this approach consists within finding the two scattering coefficients. Several analytical results exist:

• In the inviscid and adiabatic case, the results are given by Gor’kov in Ref. 3.
• Results including viscosity are given by Settnes and Bruus in Ref. 4.
• In Ref. 5, a full thermoviscous theory is given by Karlsen and Bruus. Here, the scattering is also solved off a small droplet, and the thermoviscous effects are included inside the particle as well.
• Many results in this field can be found in so-called suspension acoustics (see Ref. 6).

Equation-Based Modeling in COMSOL Multiphysics

As a first model system, we decided to determine the acoustic radiation force acting on a small spherical particle placed in an ultrasonic standing wave. As mentioned above, we opted to use the detailed perturbation approach to find the radiation force, as given in Eq. (1). The model included thermal and viscous losses in fluids and solids.

The first step was to determine the acoustic fields (the first-order field); in other words, solve the equations in the Thermoacoustics, Frequency Domain physics interface in fluids, coupled to the Thermoelasticity interface in solids. The coupling required is not predefined in COMSOL Multiphysics, but it can easily be set up in the user interface (UI). COMSOL Multiphysics will then solve the problem fully coupled.

Once the first-order fields have been established, the equations for time averages of the second-order equations can be solved. So, from here on we shall omit the \langle \cdot \rangle notation for the second-order quantities and whenever they appear in the equation, they will be understood in their time-averaged sense, \langle a_2 \rangle ­\equiv a_2.

The structure of these equations is identical to the laminar flow equations in the CFD Module, with the addition of source terms. For numerical reasons, we chose not to use the predefined equations, but rather implement our own version cast into the conservative form. This step is important for reducing the numerical noise propagating from the first-order field to the second-order solution. In the nonconservative formulation, gradients of the first-order fields need to be evaluated explicitly, thus introducing numerical noise. The equations for mass, momentum, and energy conservation are:

(4)

where it is assumed that the second-order fields are time-averaged values. Notice that the first-order fields play the role of source terms in the second-order equations.

Setting up the equations and boundary conditions was only possible using the general weak form available in COMSOL Multiphysics. The flexibility to add and define new equations was paramount to solving this problem.

The Results

We’ve covered the theory and equations. Now we can move on to looking at some nice, colorful results.

As a system example, let’s have a look at a small spherical polystyrene particle of radius from a=1 \: \mu \textrm{m} to a=40 \: \mu \textrm{m}. The particle is surrounded by air and is located in a standing 2 MHz ultrasound field, which gives \delta_\textrm{v} = 1.6 \: \mu\textrm{m} and \delta_\textrm{t} = 1.9 \: \mu\textrm{m}. Thus, for the smaller particles, boundary layer effects are expected to be significant and the lossless theory inaccurate. Moreover, for the large particles, we are outside the Rayleigh limit, k\cdot a \ll 1, with 0.04 < k\cdot a < 1.5 for the chosen particle-size range. This means that we are outside the range of the analytical theories.

The oscillating acoustic pressure field.

The animation above depicts the total harmonically oscillating acoustic pressure field p_1. The black circles represent different integration surfaces used to evaluate Eq. (1), but they are also used to control the computational mesh. The scattered field off the particle is illustrated below. This harmonic pressure oscillation leads to a net time-averaged radiation force because of nonlinear effects.

The scattered field off the particle.

The series of animations below illustrate results from a range of simulation analyses. In the top row, the figure on the left depicts temperature fluctuations T_1 around the particle. The extent of the temperature variations indicates the thermal boundary layer or penetration depth into the fluid. The figure on the right shows the deformation of the small polystyrene particle as well as the first-order temperature variations in the solid. Here, the thermal penetration depth inside the solid is again captured. Moving down to the bottom row, the figure on the left represents the axial velocity component of the scattered acoustic field, while the figure on the right indicates the radial velocity component of this field.

Results from a range of simulation analyses. Top left: Temperature variations. Top right: Deformation of the particle. Bottom left: Axial velocity component of the scattered acoustic field. Bottom right: Radial velocity component.

We can plug the first-order fields into the second-order equations — see Eq. (4) — and then solve for a local steady-state streaming flow around the particle. The pressure and velocity magnitude are shown in the next two figures, respectively. Take note of the very small scale (indicated by the color bars) for the pressure (unit: Pa) and the velocity (unit: m/s). Take, for example, the pressure relative to the incident driving wave, which has an amplitude of 3\cdot 10^{-3} \: \textrm{Pa}.

Pressure field.

Velocity field.

Finally, once the first-order and second-order fields are known, we can derive the force given in (1). In the first graph, the force is depicted as a function of \delta_\textrm{v}/a (viscous penetration depth over particle size). In the second graph, it is depicted as a function of k\cdot a. The radiation force is nondimensionalized by the corresponding force evaluated in the lossless case and in the Rayleigh limit (denoted by Gor’kov).

These graphs include three curves. The blue dashed curve is the analytical expression given by Karlsen and Bruus in the Rayleigh limit. The green dashed curve indicates the force evaluated with the lossless model (disregarding the boundary layers) given in Eq. (2). Lastly, the red curve represents our full numerical model, including thermoviscous losses, which is valid for both large and small particles.

The first graph (shown above) illustrates how the full model with viscous and thermal losses connects the lossless (\delta_\textrm{v}/a small) to the blue curve, which includes losses but is only valid for small particles (large \delta_\textrm{v}/a for constant material parameters). The second graph, which is depicted below, shows the same characteristic. The detailed thermoviscous model that we implemented connects the large particle limit (green dashed curve lossless) and the small particle limit (blue dashed curve, small k\cdot a in the Rayleigh limit).

Concluding Remarks

Attending ICSV22 gave me the opportunity to present a small portion of the research that I have worked on with Henrik Bruus and his research group at the Technical University of Denmark. In parallel with experiments and numerical modeling, the research group’s work has facilitated a deep theoretical understanding of acoustophoretic phenomena.

My presentation focused on a new numerical full thermoviscous analysis of the radiation force acting on small spherical particles (ellipsoidal particles are planned for future studies). The model agrees well with existing results in their limit of applicability (the Rayleigh limit and/or lossless limit). As this method is very general, it can be readily extended to particles of an arbitrary shape. The model can also be further expanded to include more detailed material dependencies like the fluctuation of viscosity with temperature.

Further Reading

1. “Theory and Simulation of the Acoustic Radiation Force on a Single Microparticle in an Ultrasonic Standing Wave Including Thermoviscous Effects”, M. J. H. Jensen, J. T. Karlsen, and H. Bruus.
2. “Efficient finite element modeling of radiation forces on elastic particles of arbitrary size and geometry”, P. Glynne-Jones, P. P. Mishra, R. J. Boltryk, and M. Hill.
3. “On the forces acting on a small particle in an acoustical field in an ideal fluid”, L. P. Gor’kov.
4. “Forces acting on a small particle in an acoustical field in a viscous fluid”, M. Settnes and H. Bruus.
5. “Forces acting on a small particle in an acoustical field in a thermoviscous fluid”, (accepted) J. T. Karlsen and H. Bruus.
6. Suspension Acoustics: An Introduction to the Physics of Suspensions, S. Temkin.