Numerical Simulation of Acoustic Resonance and Streaming in Droplet Acoustofluidics

TL;DR

This study uses numerical simulations to investigate how acoustic resonance weakens and streaming begins in droplet acoustofluidics, revealing the influence of phase properties on energy and streaming phenomena.

Contribution

It introduces a numerical approach to analyze resonance weakening and streaming onset in droplet acoustofluidics based on phase property variations.

Findings

Resonance weakening observed with increased wall oscillation amplitude.

Streaming onset depends on relative acoustic properties of phases.

Implications for acoustic focusing inside droplets.

Abstract

This paper reports the phenomenon of resonance weakening and streaming onset in two phase acoustofluidics by performing numerical simulations of a capillary droplet suspended in a microfluidic chamber. The simulations show that depending on the relative acoustic properties of the two phases, it is possible to observe (i)~the decrease in the total acoustic energy as the oscillation amplitude at the wall increases, and (ii)~the onset of acoustic streaming. The impact of these findings in terms of acoustic focusing inside droplets is also discussed.

Tables

Table 1. Table 1: List of the physical properties used in the numerical simulations.

	${\rho }_0\left[\mathrm{kg}/m^3\right]$	$c_0\left[m/s\right]$	${\eta }_0\left[\mathrm{mPa}\cdot s\right]$	${\gamma }_{W,x}\left[\mathrm{mN}/m\right]$
Water	$997$	$1497$	$1.0$	$-$
Olive Oil	$910$	$1450$	$61.0$	$16.0$
HFE-7500	$1619$	$660$	$1.2$	$41.0$
Silicon Oil	$950$	$1100$	$19.0$	$32.0$

Equations

$${\partial }_t\rho$$

$${\partial }_t\left(\rho \mathbf{v}\right)$$

$$\mathbf{\Pi }=\rho \mathbf{v}\mathbf{v}-\mathbf{\sigma },$$

$$\mathbf{\sigma }=-\left(p-\lambda \nabla \cdot \mathbf{v}\right)\mathbf{I}+\eta \left[\nabla \mathbf{v}+{\left(\nabla \mathbf{v}\right)}^T\right],$$

$${\mathbf{v}}_{\mathrm{w}}(\pm W_c/2,z)=[U(z),0{]}^T{e}^{-i\omega t},$$

$$\mathbf{v}(y,\pm H_c/2)=0,$$

$$0$$

$${\mathbf{f}}_{\mathrm{st}}=\gamma \left({\nabla }_{\mathrm{s}}\cdot \hat{\mathbf{n}}\right)\hat{\mathbf{n}}.$$

$${\nabla }_{\mathrm{s}}=\left(\mathbf{I}-\hat{\mathbf{n}}\hat{\mathbf{n}}\right)\cdot \nabla ,$$

$${\mathcal{L}}_{\mathrm{hd}}^{\mathrm{int}}=-\gamma \mathrm{Tr}\left[\left(\mathbf{I}-\hat{\mathbf{n}}\hat{\mathbf{n}}\right)\cdot \nabla {\mathbf{v}}_0\right]+{\mathbf{f}}_{\mathrm{ac}}\cdot {\mathbf{v}}_0,$$

$$0$$

$${\partial }_t\left({\rho }_0{\mathbf{v}}_0\right)$$

$${\mathbf{\Pi }}_0={\rho }_0{\mathbf{v}}_0{\mathbf{v}}_0-{\mathbf{\sigma }}_0,$$

$${\mathbf{\sigma }}_0=-p_0\mathbf{I}+{\eta }_0\left[\nabla {\mathbf{v}}_0+{\left(\nabla {\mathbf{v}}_0\right)}^T\right].$$

$${\mathbf{v}}_0=0.$$

$$0$$

$${\mathbf{f}}_{\mathrm{ac}}=-[[⟨{\mathbf{\Pi }}_2⟩]]\cdot \hat{\mathbf{n}}.$$

$$⟨{\mathbf{\Pi }}_2⟩=\frac{⟨\nabla p_1\nabla p_1⟩}{{\rho }_0{\omega }^2}+\left(\frac{⟨p_1^2⟩}{2{\rho }_0{c}_0^2}-\frac{⟨⟩\nabla p_1{}^2}{2{\rho }_0{\omega }^2}\right)\mathbf{I}$$

$$⟨{\mathbf{\Pi }}_2⟩=\frac{\nabla {\hat{p}}_1\nabla {\hat{p}}_1}{2{\rho }_0{\omega }^2}+\left(\frac{{\hat{p}}_1^{\ast }{\hat{p}}_1}{4{\rho }_0{c}_0^2}-\frac{\nabla {\hat{p}}_1^{\ast }\cdot \nabla {\hat{p}}_1}{4{\rho }_0{\omega }^2}\right)\mathbf{I}$$

$$\frac{{\hat{p}}_1}{{\rho }_0{c}_0^2}+\nabla \cdot \left(\frac{\nabla {\hat{p}}_1}{{\rho }_0{\omega }^2}\right)=0,$$

$$\nabla {\hat{p}}_1\cdot \hat{\mathbf{n}}=\mathrm{i}{\rho }_0\omega [U(z),0{]}^T,$$

$$\nabla {\hat{p}}_1\cdot \hat{\mathbf{n}}=0,$$

$$[[{\hat{p}}_1]]$$

$$[[{\rho }_0^{-1}\nabla {\hat{p}}_1]]\cdot \hat{\mathbf{n}}$$

$${\mathcal{L}}_{\mathrm{ac}}(\mathbf{x},t)=\frac{\nabla {\hat{p}}_1^{\ast }\cdot \nabla {\hat{p}}_1}{4{\rho }_0{\omega }^2}-\frac{{\hat{p}}_1^{\ast }{\hat{p}}_1}{4{\rho }_0{c}_0^2},$$

$${\mathcal{P}}_{\mathrm{ac}}=\frac{1}{2}{\hat{p}}_1^{\ast }{\hat{\mathbf{v}}}_1\cdot \hat{\mathbf{n}},$$

$${\mathcal{E}}_{\mathrm{ac}}(\mathbf{x},t)=\frac{\nabla {\hat{p}}_1^{\ast }\cdot \nabla {\hat{p}}_1}{4{\rho }_0{\omega }^2}+\frac{{\hat{p}}_1^{\ast }{\hat{p}}_1}{4{\rho }_0{c}_0^2},$$

$$E_{\mathrm{ac}}(t)={\int }_V{\mathcal{E}}_{\mathrm{ac}}(\mathbf{x},t)\mathrm{d}V,$$

$${\epsilon }_{\mathrm{ac}}(t)=\frac{E_{\mathrm{ac}}(t)}{E_{\mathrm{ac}}(0)},$$

$$V_{\mathrm{droplet}}={\int }_{\mathrm{droplet}}\mathrm{d}V$$

$${\mathbf{x}}_{\mathrm{droplet}}=\frac{1}{V_{\mathrm{droplet}}}{\int }_{\mathrm{droplet}}\mathbf{x}\mathrm{d}V$$

Full text

I Introduction

Acoustofluidics is an established microfluidics technology that combining ultrasounds and a carrier flow allows for a fine manipulation of particle at microscale Evander and Nilsson (2012); Augustsson and Laurell (2012); Nordin and Laurell (2012); Manneberg et al. (2008); Augustsson et al. (2012); Ding et al. (2014); Petersson et al. (2007).

Recently, some researchers attempted to combine acoustofluidics and droplet acoustofluidics for particle focusing inside droplets. The results for standard fluids used in droplet microfluidics, i.e. fluorinated oil HFE-7500, have been subjected of controversies. Indeed, it seems that other organic phases have shown better performance for particle focusing inside the droplets.

The present paper employs numerical simulations to give a first insight on the phenomena involved in acoustic resonance and streaming in droplet acoustofluidics. Specifically, simulations of a water droplet in HFE-7500 and olive oil have been performed.

II Model Setting

The geometry used in the simulation is that reported in Fig. 1. A water droplet (suffix “’w’ means water properties) is suspended in an organic phase (suffixes “OO” and “HFE” refer to olive oil and HFE respectively) within a microfluidic chamber. The vertical channel walls oscillate with a predetermined displacement amplitude.

Oscillation Amplitude at Chamber Walls.

Since the piezo and the material embedding the microchannel are neither modeled nor their response is simulated, there is the need to assume the acoustic pressure amplitude at the resonance frequency. This implies also that there is no feedback to the transducer, namely no change in the overall impedance. The driving pressure rescaling is performed by computing the resonance acoustic frequency and then rescaling the driving amplitude at the boundary proportionally to the desired value of the acoustic pressure amplitude $p_{\mathrm{ac}}$, defined as the maximum of the acoustic pressure. The hypothetical acoustic pressure amplitudes considered in the simulations are in the range $p_{\mathrm{ac}}=0.01-5\,\mathrm{MPa}$.

II.1 Governing Equations

The governing equations are the mass density conservation and the linear momentum density balance,

[TABLE]

where $\rho$ is the fluid density, $\bm{v}$ is the fluid velocity, $\bm{f}_{\mathrm{v}}$ is a volume force (there is a plan to include gravity, computing the initial configuration for the droplet and then run the time-dependent simulations), and

[TABLE]

is the momentum (density) flux. In the latter relation, $\rho\bm{v}\bm{v}$ is the advective momentum flux while $\bm{\sigma}$ is the stress tensor

[TABLE]

in which $p$ is the pressure, $\eta$ is the shear viscosity, and $\lambda$ is the dilational viscosity (values for this parameter are difficult to find, we consider either $\lambda=0$ or by Stokes’ approximation $\lambda\simeq\tfrac{2}{3}\eta$).

Equations need to be equipped with the boundary conditions. Here a rectangular closed domain (see Fig. 1) with two vertical walls at $y=\pm W_{c}/2$ is considerd. The vertical walls oscillate at angular frequency $\omega$, and thus it can write for the wall velocity $\bm{v}_{\mathrm{w}}$

[TABLE]

where $U(z)$ is chosen such that the simulation yields the desired pressure amplitude, that is $p_{\mathrm{ac}}$. At the bottom and top walls the rigid wall conditions

[TABLE]

are assigned.

II.2 Interface Conditions

At the interface between the continous phase and disperse phase (the droplet) the interface conditions

[TABLE]

are considered ($\bm{\hat{n}}$ is the outward-pointing normal). The surface tension force is given by

[TABLE]

The symbol $\bm{\nabla}_{\mathrm{s}}$ denotes the surface nabla operator,

[TABLE]

where $\bm{I}$ is the unit tensor ($[\bm{I}]_{ij}=\delta_{ij}$).

Weak form for Surface Tension at the Interface.

In the presence of surface tension force and acoustic force, the weak form contribution at the interface is given by

[TABLE]

which is true just for a closed interface. Note that there is no need to compute the interface curvature since in the weak form this will not appear.

II.3 Time-scale Splitting

It is assumed that the hydrodynamic time-scale and the acoustic time-scale are well separated and correspond to a slow time-scale and a fast time-scale, respectively. The slow-time scale dynamics is described by means of the incompressibility constraint and the linear momentum density balance for the time-average density $\rho_{0}$ and velocity $\bm{v}_{0}$,

[TABLE]

The time-averaged momentum flux $\bm{\Pi}_{0}$ is

[TABLE]

where

[TABLE]

For the slow time-scale we consider the no-slip conditions at the channel walls

[TABLE]

The slow time-scale interface conditions are

[TABLE]

where $\bm{f}_{\mathrm{ac}}$ is the acoustic force at the interface

[TABLE]

The time-average of the second-order momentum density flux is given by

[TABLE]

or in terms of complex-valued pressure amplitude $\hat{p}_{1}$

[TABLE]

where the amplitude of the first-order acoustic pressure, i.e. $p_{1}(\bm{x},t)=\hat{p}_{1}(\bm{x},\omega)\,e^{-i\,\omega\,t}$, is the solution of the Helmholtz equation

[TABLE]

with the boundary condition

[TABLE]

at the vertical walls and

[TABLE]

at the bottom and top walls. The fast-time scale interface conditions are

[TABLE]

III Results

Table LABEL:tab:physprop reports the physical properties used in the num erical simulations.

Figure 2 shows the resonance spectrum (a) in the case of a water droplet suspended in olive oil, silicone oil, and HFE oil (see caption for color coding), by plotting the acoustic energy $E_{\mathrm{ac}}(0)$ as function of the frequency $f$. The three resonance frequencies set in the beginning of the time-dependent simulations were (put the resonance frequencies here).

Figure 2 shows the ratio between the acoustic energy when the droplet has reach the final stage and the initial acoustic energy, that is $\varepsilon^{\infty}_{\mathrm{ac}}$, as function of the rescaled driving acoustic pressure at the walls $p_{\mathrm{ac}}$. It has been found that as the driving pressure increases the final acoustic energy is less than that expected for a immutable geometric configuration, in the case of olive oil and silicone oil. In the case of HFE the acoustic energy has shown a small increase before to reach a value for the driving pressure $p_{\mathrm{ac}}$ beyond that the droplet dynamics is unstable.

In the case of the stable simulations for HFE/water system, it has been found that there is an onset of acoustic streaming.

The evidence of this is reported in figure 3 that shows the ratio between the kinetic energy and the final kinetic energy as function of the time (a), together with the average fluid velocity at the final time.

The phenomenon of acoustic streaming onset has not been observed for OO/water and SO/water system.

The evidence of this fact has been reported in figure 4 where the together with the rescaled pressure field (color shading), also the normalized acoustic streaming velocity field is reported (blue arrows).

IV Discussion

As one studies the behavior of the acoustic energy for a closed system where a droplet of water is immersed in an organic phase, two different different phenomenologies can be observed. These depend on on the relative acoustic properties of the two phases.

IV.1 Acoustic Energy

For OO/water systems, as the driving acoustic pressure increases the energy appears to be initially constant as the interface tension force prevents the particle deformation. However, when the acoustic force at the interface starts to increase and competing with the interface tension force the droplet is deformed. As this deformation appears, the driving frequency does not match the resonance frequency of the new geometric configuration and the acoustic energy descreases.

For HFE/water systems, the observations suggest that as the driving pressure amplitude increases, the energy increases until an unstable state (most likely corresponding to a new resonance frequency for the deformed system) is reached. The simulations could not be performed beyond a certain value of the acoustic driving pressure.

IV.2 Acoustic Streaming

For OO/Water systems, the relative physical properties of the two fluids have shown a suppression of the acoustic streaming.

For HFE/Water systems, it has been observed that there is an onset of the acoustic streaming. This fact deserve more investigations.

V Conclusions

The impact of having a capillary droplet enclosed in a microfluidic and acoustically driven microchamber has been studied.

Two behaviors has been observed and these depend on the relative acoustic properties of the disperse phase (the droplet) and the continous phase (the surrounding fluid).

In the case of a water droplet suspended in olive oil, it has been observed both a weakening of the resonance and streaming suppression. In the case of a water droplet suspended in HFE oil, it has been observed an increase in the acoustic energy and the onset of streaming.

Further investigations are considered since it is important to provide a meaningful explanation of the observed behaviors and how these impact on particle/cell manipulation inside droplets.

Appendix A Indicators

(short scale time-averaged)Acoustic Lagrangian Density

[TABLE]

Acoustic Power-Input Density

[TABLE]

Acoustic Energy Density (as indicator of acoustophoresis)

[TABLE]

the total acoustic energy is then

[TABLE]

as the initial configuration is that used to compute the resonance frequency we can define the ratio between the final acoustic energy and the initial acoustic energy, that depends on $p_{\mathrm{ac}}$

[TABLE]

note that for an “immutable” geometric configuration this ratio must equals $1$ for any $t$. We define $\varepsilon_{\mathrm{ac}}^{\infty}=E_{\mathrm{ac}}(\infty)/E_{\mathrm{ac}}(0)$.

A.1 Indicators

A.1.1 Assessment of good simulations

• •

simulation of small droplets dislocation towards to/far from the pressure node (depending on the contrast factor)

• •

comparison with standard theory of particle motion in acoustic field (give a try to the shift in resonance when the particle/droplet undergoes focusing, it can be interested to show this)

Analysis of Droplet Motion

Droplet Volume (in 2D this measures an area)

[TABLE]

this in reliable simulations must be a constant, or at least we admit $\pm 5\%$ in variation (or less depending on the computational efforts). The averaged droplet position can be computed by considering the integral

[TABLE]

when we validate the simulations, we compare this quantity against the theoretical droplet trajectory given by integration of the equation of motion for a spherical particle subjected to acoustic radiation forces. In the case of the centered droplet we admit a variation pf $\pm 5\%$ (maybe less) from the centered position. Acoustic force on the droplet

[TABLE]

to be compared with the known expressions in the case of “particle-like” droplet.

A.1.2 Phenomenon Analysis

A.1.3 Acoustics

A.1.4 Hydrodynamics

Hydrodynamic Lagrangian Density

[TABLE]

rate of change in the hydrodynamic kinetic energy density

[TABLE]

Hydrodynamic kinetic energy density,

[TABLE]

needed to evaluate the presence/absence of a velocity field in the model.

(Rate of) viscous dissipation density

[TABLE]

Total kinetic energy density

[TABLE]

[TABLE]

characteristic system velocity

[TABLE]

where $m=\int_{V}\rho(\bm{x},t)\mathrm{d}V=\mathrm{const}$. We also define $V^{\infty}=V(t\rightarrow\infty)$. Hydrodynamic power input density (at the interface between the two fluids)

[TABLE]

Lagrangian density for the time-dependent model

Note that

[TABLE]