2.1.1 Periodicity and the unit cell

We mostly consider flows in 3-D periodic systems. To this end, we introduce a triclinic unit cell $\unicode[STIX]{x1D6E4}\subset \mathbb{R}^{3}$ that is spanned by the basis $\{\boldsymbol{a}^{(1)},\boldsymbol{a}^{(2)},\boldsymbol{a}^{(3)}\}$ as shown in figure 1(a). The three basis vectors $\boldsymbol{a}^{(i)}\in \mathbb{R}^{3}$ , $i=1,2,3$ , define a right-handed coordinate system. In the most general case they are not required to be of unit length or orthogonal and might depend on time, although the latter will not be explicitly considered here. Nevertheless, even static but skewed bases can be useful in practice, for example to model so-called deterministic lateral displacement devices (Huang et al. Reference Huang, Cox, Austin and Sturm2004; Krüger, Holmes & Coveney Reference Krüger, Holmes and Coveney2014; Zhang et al. Reference Zhang, Henry, Gompper and Fedosov2015). The surface $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ of the unit cell is not included in the open set  $\unicode[STIX]{x1D6E4}$ . We will denote the volume by $V_{\unicode[STIX]{x1D6E4}}$ .

Figure 1.

(a) Example of a triclinic unit cell $\unicode[STIX]{x1D6E4}$ (thick lines), together with 11 replicas (thin grey lines). In this example, the two basis vectors $\boldsymbol{a}^{(1)}$ and $\boldsymbol{a}^{(2)}$ lie in the $x$ – $y$ -plane, but $\boldsymbol{a}^{(2)}$ is not parallel to the $y$ -axis. (b) Two-dimensional sketch of the general three-dimensional problem. $\unicode[STIX]{x1D6E4}\subset \mathbb{R}^{3}$ contains the whole unit cell (everything within the dashed border $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ ).  $\unicode[STIX]{x1D6FA}\subset \unicode[STIX]{x1D6E4}$ is the ambient fluid with dynamic viscosity  $\unicode[STIX]{x1D707}$ (light blue). As an example, two capsule-like entities ( ${\mathcal{C}}_{1}$ , ${\mathcal{C}}_{2}$ ), two volume-changing bubbles ( ${\mathcal{B}}_{1}$ , ${\mathcal{B}}_{2}$ ) and a wall ${\mathcal{W}}_{1}$ are shown. The normal vectors always point into $\unicode[STIX]{x1D6FA}$ . The arrows on the left represent an imposed flow.

To cast the concept of periodicity into mathematical terms, we introduce by

some grid vector with some grid index $\unicode[STIX]{x1D736}\in \mathbb{Z}^{3}$ . The unit cell $\unicode[STIX]{x1D6E4}$ is assumed to be replicated infinitely throughout space by displacing it with all possible grid vectors. We will show in appendix B that actually any of these boxes with any origin can be chosen as $\unicode[STIX]{x1D6E4}$ . Furthermore, a general function $f:\mathbb{R}^{3}\rightarrow \mathbb{R}$ is said to periodic if it satisfies

The reciprocal (Fourier) unit cell is spanned by the reciprocal basis vectors $\boldsymbol{b}^{(j)}\in \mathbb{R}^{3}$ , $j=1,2,3$ , defined by $\boldsymbol{a}^{(i)}\boldsymbol{\cdot }\boldsymbol{b}^{(j)}=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}_{ij}$ . Some general Fourier grid vector is then denoted by

with $\unicode[STIX]{x1D73F}\in \mathbb{Z}^{3}$ .

2.1.2 System components

The general system considered in this section is depicted in figure 1(b). We introduce the ambient fluid $\unicode[STIX]{x1D6FA}\subset \unicode[STIX]{x1D6E4}$ which represents the open set of the space within the unit cell  $\unicode[STIX]{x1D6E4}$ but outside of any immersed object. The fluid is assumed to have a constant dynamic viscosity $\unicode[STIX]{x1D707}$ and is governed by the usual Stokes equation and incompressibility condition. Furthermore, we have $N_{O}$ objects $O_{k}\subset \unicode[STIX]{x1D6E4}$ , $k=1,\ldots ,N_{O}$ in the unit cell as further detailed below. For now we assume that all of them are completely located within $\unicode[STIX]{x1D6E4}$ and relax this requirement later on (which is required for dynamic simulations of dense suspensions, see appendix B). No object shall overlap with or contain any other object. The 2-D surfaces of the objects embedded in 3-D space will be denoted by the symbol ‘ $\unicode[STIX]{x2202}$ ’, e.g.  $\unicode[STIX]{x2202}O_{k}$ . We also assign a velocity $\boldsymbol{u}(\boldsymbol{x})$ to each point $\boldsymbol{x}\in \bar{\unicode[STIX]{x1D6E4}}$ , where the bar represents the closure of the set.

Regarding the immersed objects, we consider three different types. First, deformable particles that can be used to mimic ‘capsule-like’ entities (such as vesicles, drops or red blood cells). The $i$ th capsule will be denoted by the open set ${\mathcal{C}}_{k}\subset \unicode[STIX]{x1D6E4}$ with $k=1,\ldots ,N_{{\mathcal{C}}}$ , where $N_{{\mathcal{C}}}$ are the number of capsules. Their inside is filled with some Stokesian fluid that has a dynamic viscosity of $\unicode[STIX]{x1D707}\unicode[STIX]{x1D706}_{{\mathcal{C}}_{k}}$ , where $\unicode[STIX]{x1D706}_{{\mathcal{C}}_{k}}>0$ is the viscosity ratio between the inner and outer fluids. Their volume is conserved and as such the net flux $Q_{{\mathcal{C}}_{k}}$ into or out of the capsules is zero. Deformations are governed by the jump of the traction $\unicode[STIX]{x0394}\boldsymbol{f}$ across their surfaces $\unicode[STIX]{x2202}{\mathcal{C}}_{k}$ . This traction jump is calculated from an appropriate constitutive law that determines the actual object properties. The velocity across the interfaces is assumed to be continuous.

Second, $N_{{\mathcal{W}}}$ non-closed objects ${\mathcal{W}}_{k}$ , $k=1,\ldots ,N_{{\mathcal{W}}}$ , may exist that can be used to model deformable (fixed traction jump) or rigid (prescribed velocity) walls. Note that we set ${\mathcal{W}}_{k}=\unicode[STIX]{x2202}{\mathcal{W}}_{k}$ , i.e. these objects do not have some ‘inside’. The velocity across the walls is assumed to be continuous.

Third, we introduce our novel ingredient, $N_{{\mathcal{B}}}$ bubble-like objects ${\mathcal{B}}_{k}$ , $k=1,\ldots ,N_{{\mathcal{B}}}$ whose volumes are allowed to change. Contrary to capsules their inside is filled with some compressible fluid such as air that has a much lower viscosity than the ambient fluid. Therefore, the detailed flow field inside the bubbles will not be considered and, instead of the Stokes equations, the essential model assumption for this inner fluid is simply

Here, $c_{k}$ are some prescribed constants, meaning that the particles expand or contract homogeneously over their whole interior and that any inhomogeneities such as eddies are considered to vanish instantaneously. But note that $c_{k}$ may depend on the time allowing for oscillating bubbles. As shown below (2.31), $c_{k}$ is related to the net flux $Q_{{\mathcal{B}}_{k}}$ into or out of the bubble by $c_{k}=Q_{{\mathcal{B}}_{k}}/V_{{\mathcal{B}}_{k}}$ , where $V_{{\mathcal{B}}_{k}}$ is the $k$ th bubble’s volume. Moreover, the surface deformation of the bubbles is governed by the prescribed traction at the outside of the surfaces, as further elaborated in § 2.4.

For convenience, some arbitrary object of any type will be denoted by $O_{k}$ . We further define $N_{O}:=N_{{\mathcal{C}}}+N_{{\mathcal{B}}}+N_{{\mathcal{W}}}$ as well as $O:={\mathcal{C}}\cup {\mathcal{B}}\cup {\mathcal{W}}$ where ${\mathcal{C}}$ / ${\mathcal{B}}$ / ${\mathcal{W}}$ mean the unions of all capsules/bubbles/walls. All in all, we have $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6FA}\cup \bar{O}$ . Furthermore, unit normal vectors are denoted by the symbol $\boldsymbol{n}$ . We choose the convention that the normal vector of every capsule and bubble points into $\unicode[STIX]{x1D6FA}$ . The normal vector of walls may point in any of the two possible directions. The normal vector on $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ shall point into $\unicode[STIX]{x1D6FA}$ .

2.2.1 The general boundary integral equation

In order to introduce the required notation and equations for later on, we start from the standard Stokes and continuity equations for flows at low Reynolds numbers which corresponds to the typical situation encountered in the microcirculation (Misbah Reference Misbah2012; Freund Reference Freund2014):

for $\boldsymbol{x}\in \unicode[STIX]{x1D6FA}$ . Here, $P$ is the pressure, $\unicode[STIX]{x1D6FB}^{2}$ the usual 3-D Laplace operator and $\boldsymbol{u}$ the fluid velocity. The inside of capsules is governed by analogous expressions, but with a different viscosity in general. Body forces such as gravity will be neglected, but can be easily incorporated via an effective pressure (Pozrikidis Reference Pozrikidis1992, equation (1.2.9)).

Next, we introduce the traction $\boldsymbol{f}$ acting on the surface $\unicode[STIX]{x2202}O_{k}$ of some object $O_{k}$ via

where the fluid stress tensor is defined by

Summation over repeated Cartesian components is implied throughout this work. One needs to take the limit of the stress tensor onto the surface for the evaluation of (2.7). This limit can be taken from both sides, leading to the traction on the outside (  $\boldsymbol{f}^{+}$ , limit from $\unicode[STIX]{x1D6FA}$ onto the surface) and the inside (  $\boldsymbol{f}^{-}$ , defined with the outer normal vector) of a closed object’s surface. The difference is the traction jump

which is the major quantity coupling the flow with the surface mechanics. For a wall, $\unicode[STIX]{x0394}\boldsymbol{f}$ is the difference of the tractions on its two sides.

Equations (2.6) can be efficiently and accurately solved via the boundary integral method (Pozrikidis Reference Pozrikidis1992, Reference Pozrikidis2001; Freund Reference Freund2014). Assuming that all objects are located within the unit cell $\unicode[STIX]{x1D6E4}$ (for objects crossing the unit cell boundary, see appendix B), one can derive the boundary integral (BI) equation (compare Pozrikidis Reference Pozrikidis1992, chap. 2.3, pp. 37 and 143)

Note that the evaluation point $\boldsymbol{x}_{0}$ is in the ambient fluid and not within any object or on any surface. The single-layer integral (or single-layer potential) over some surface  $S$ is defined as

and the double-layer integral is

These integrals exist for $\boldsymbol{x}_{0}\in \mathbb{R}^{3}$ , notably in the improper sense if $\boldsymbol{x}_{0}\in S$ or if $\boldsymbol{x}_{0}$ is on some periodic image of $S$ (Kress Reference Kress2014, chap. 6). The improper double-layer integral is in this case sometimes also denoted as a principal value integral (Pozrikidis Reference Pozrikidis1992, p. 27).

$\boldsymbol{F}$ is a ‘unified traction’ which represents either the traction $\boldsymbol{f}^{+}$  on the outside in case of the bubbles, or the traction jump $\unicode[STIX]{x0394}\boldsymbol{f}$ in the case of the capsules and the walls. For notational convenience we have abbreviated $\boldsymbol{f}:=\boldsymbol{f}^{+}$ in the last line of (2.10). Moreover, an ‘effective’ viscosity ratio is defined as

Thus, if , the corresponding terms in the second sum are always absent (regardless of the value of ${\mathcal{K}}_{\unicode[STIX]{x2202}O_{q}}\boldsymbol{u}$ ).

The wall contributions in (2.10) can be derived by considering as a starting point an object with a finite thickness. Taking the limit as the thickness tends to zero, integrals that previously went over the distinct opposing sides of the wall now go over essentially the same surface, except that the integrands are still evaluated on their respective sides. The normal vectors on these two sides are perfectly antiparallel. Additionally using the continuity of the velocity, one finds that the double-layer integral drops out (represented by ). Furthermore, for the single-layer integral we have $\int _{\unicode[STIX]{x2202}{\mathcal{W}}^{+}\cup \unicode[STIX]{x2202}{\mathcal{W}}^{-}}f_{i}^{+}G_{ij}\,\text{d}S=\int _{\unicode[STIX]{x2202}{\mathcal{W}}^{+}}(f_{i}^{+}-f_{i}^{-})G_{ij}\,\text{d}S\equiv \int _{\unicode[STIX]{x2202}{\mathcal{W}}}\unicode[STIX]{x0394}f_{i}G_{ij}\,\text{d}S$ where $\unicode[STIX]{x2202}{\mathcal{W}}^{\pm }$ denotes the two sides of the zero-thickness wall, and (2.9) was used to introduce the traction jump (Pozrikidis Reference Pozrikidis1992, p. 37). The integration $\int _{\unicode[STIX]{x2202}{\mathcal{W}}}$ goes over only one side of the wall. Walls spanning through the whole domain $\unicode[STIX]{x1D6E4}$ can be introduced by a similar limiting procedure.

Note that the BI equation (2.10) is actually valid for any Green’s functions, no matter if periodic or not. However, in order to be able to compute anything in practice, they should follow the general symmetries of the system. As we want to implement a periodic domain, we will next introduce the appropriate expressions.

2.2.2 Appropriate Green’s functions for a 3-D periodic domain

Green’s functions $G_{ij}$ and $p_{j}$ for the velocity and pressure, respectively, are obtained by solving the singularly forced Stokes equation. It is possible to derive specialized Green’s functions such as for singly (Pozrikidis Reference Pozrikidis1996) or doubly (Pozrikidis Reference Pozrikidis1996; Greengard & Kropinski Reference Greengard and Kropinski2004; Cortez & Hoffmann Reference Cortez and Hoffmann2014) periodic domains or with incorporated stationary walls (Blake Reference Blake1971; Liron & Mochon Reference Liron and Mochon1976; Staben, Zinchenko & Davis Reference Staben, Zinchenko and Davis2003; Janssen & Anderson Reference Janssen and Anderson2008). However, we want to be able to apply our method for general channel geometries, and therefore employ the standard Green’s functions for a 3-D periodic domain. According to Hasimoto (Reference Hasimoto1959) and Pozrikidis (Reference Pozrikidis1996), they are given by

for the velocity ( $G_{ij}$ ) and the pressure ( $p_{j}$ ), respectively. $\boldsymbol{x}_{0}\in \bar{\unicode[STIX]{x1D6E4}}$ is arbitrary but fixed. The sums go over all possible Fourier grid vectors $\boldsymbol{k}^{(\unicode[STIX]{x1D73F})}$ as defined in (2.4) with $\unicode[STIX]{x1D73F}\in \mathbb{Z}^{3}\setminus \{\mathbf{0}\}$ , and $\text{i}$ is the imaginary unit. $G_{ij}$ satisfies the incompressibility condition,

Furthermore, $\langle G_{ij}\rangle _{\unicode[STIX]{x1D6E4}}:=(1/V_{\unicode[STIX]{x1D6E4}})\int _{\unicode[STIX]{x1D6E4}}G_{ij}(\boldsymbol{x},\boldsymbol{x}_{0})\,\text{d}x^{3}$ is independent of $\boldsymbol{x}_{0}$ and constitutes a free parameter representing an imposed average net flow for the Green’s function system. We choose the coordinate system such that it is zero (Zinchenko & Davis Reference Zinchenko and Davis2000):

Note that the full system can nevertheless have an imposed flow which will be introduced in § 2.2.4 via double-layer integrals over the unit cell boundary.

Combining both of the above equations (via Pozrikidis Reference Pozrikidis1992, equation (2.1.8)), the stresslet is found to be

$T_{ijl}$ contains a linear (non-periodic) part, whereas the periodic part is

The periodicity holds for both arguments, as well as for the Stokeslet:

This also implies that $\boldsymbol{x}_{0}$ can actually be located within almost the whole space $\mathbb{R}^{3}\setminus \{\boldsymbol{x}\}$ .

Equations (2.14) and (2.18) are the Green’s functions used below. The above given forms, however, are not computable in practice due to their slow convergence. Dramatic speedups are achieved by using the Ewald decomposition (Hasimoto Reference Hasimoto1959). The final formulas are given by Zhao et al. (Reference Zhao, Isfahani, Olson and Freund2010).

We will additionally need two general relations. Hence, let $O$ be some object with a closed surface $\unicode[STIX]{x2202}O$ , and let $\{O^{(\unicode[STIX]{x1D736})}\}$ be the set of $O$ that is offset with all possible periodic grid vectors $\boldsymbol{X}^{(\unicode[STIX]{x1D736})}$ from (2.2). Then, for some proper Green’s functions:

2.2.3 Boundary conditions for a 3-D periodic domain

As we wish to simulate a 3-D periodic domain with the BI equation (2.10), we stipulate as a boundary condition that the velocity $\boldsymbol{u}$ shall be periodic, i.e. 

Note that in principle $\boldsymbol{u}$ may contain a linear component, as long as the arrangement of the objects and their images retain a triclinic periodicity. This could be used to simulate an overall linear shear flow (Phan-Thien, Tran-Cong & Graham Reference Phan-Thien, Tran-Cong and Graham1991; Pozrikidis Reference Pozrikidis1993, Reference Pozrikidis1999; Li, Zhou & Pozrikidis Reference Li, Zhou and Pozrikidis1995; Loewenberg & Hinch Reference Loewenberg and Hinch1996). However, this also usually implies that the basis vectors $\boldsymbol{a}^{(i)}$ of the unit cell have to be time-dependent, leading to additional problems for longer simulation times (Blawzdziewicz Reference Blawzdziewicz, Andrea and Gretar2007, p. 221).

As for the pressure Green’s function from (2.15), it is possible to have a non-zero pressure gradient over the unit cell that drives a certain average flow $\langle \boldsymbol{u}\rangle _{\unicode[STIX]{x1D6E4}}$ . Hence, we can decompose the pressure as $P(\boldsymbol{x})=\langle \unicode[STIX]{x1D735}P\rangle _{\unicode[STIX]{x1D6E4}}\boldsymbol{\cdot }\boldsymbol{x}+\breve{P}(\boldsymbol{x})$ for $\boldsymbol{x}\in \mathbb{R}^{3}$ where $\breve{P}$ denotes the periodic part, and

is the average over the whole unit cell. Using (2.8) and (2.7), we find for the traction (also compare Zick & Homsy Reference Zick and Homsy1982, equation (2.23))

where $\breve{\boldsymbol{f}}$ denotes the periodic part.

We finally remark that the BI equation (2.10) can be used with any Green’s function to simulate a periodic domain, provided that the periodicity condition (2.23) is enforced. Using the periodic Green’s functions from (2.14) and (2.18) is merely a convenient choice since the integrals over $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ then yield simple expressions, as shown next.

2.2.4 Computing the integrals over the unit cell’s surface for a periodic domain

Single-layer integral. With the help of the decomposition (2.25), the periodicity of the involved quantities (due to (2.20a ) and (2.23)) and the divergence theorem together with (2.17) an (2.16), one can show (Zick & Homsy Reference Zick and Homsy1982)

Double-layer integral. It is tempting to assume that the double-layer integral ${\mathcal{K}}_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}}\boldsymbol{u}$ over the unit cell surface $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ is also zero. This, however, will turn out to be wrong if a net flow is imposed and/or if volume-changing objects are included. The latter will lead to an additional novel contribution to the equations. Following Zick & Homsy (Reference Zick and Homsy1982), the first step for the evaluation is to use the decomposition (2.18). The integral term containing $\breve{T}_{ijl}$ vanishes due to periodicity (compare (2.23) and (2.20b ), and note that normal vectors on opposing sides of $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ are antiparallel). The other term is treated by adding and subtracting integrals over the objects, resulting in

Integrals over walls give zero contributions since the velocity $\boldsymbol{u}$ is continuous across their surface, but the normal vectors on opposite sides have different signs. Using the continuity of the normal velocity across the interfaces and applying the divergence theorem while watching out for the normal vector convention (always into the ambient fluid $\unicode[STIX]{x1D6FA}$ ), one obtains

The integrals in the first line can be combined to $\int _{\unicode[STIX]{x1D6E4}}\boldsymbol{u}\,\text{d}x^{3}=V_{\unicode[STIX]{x1D6E4}}\langle \boldsymbol{u}\rangle _{\unicode[STIX]{x1D6E4}}$ , with the average defined in (2.24) (walls are nullsets). Moreover, the divergence of the velocity vanishes in $\unicode[STIX]{x1D6FA}$ and the capsules ${\mathcal{C}}$ because of (2.6b ) (Stokesian fluids). Furthermore, the last term is absent in existing formulations without volume-changing objects, but here it is non-zero in general and therefore requires special attention.

A more usable form of this last term may be obtained by using the model from (2.5). For a particular bubble ${\mathcal{B}}_{k}$ , $k=1,\ldots ,N_{{\mathcal{B}}}$ , we immediately find

where we have defined the geometric centroid

$V_{{\mathcal{B}}_{k}}$ is the bubble volume. We obtain a connection between $c_{k}$ to the more intuitive flux  $Q_{{\mathcal{B}}_{k}}$ out of or into a bubble by computing

Putting it all together, the double-layer integral over $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ hence becomes

Similar to $\langle G_{ij}\rangle _{\unicode[STIX]{x1D6E4}}$ for the Green’s function from (2.14), the average velocity $\langle u_{j}\rangle _{\unicode[STIX]{x1D6E4}}$ is a free parameter that can be used to drive a flow through the system (Zhao et al. Reference Zhao, Isfahani, Olson and Freund2010). Since the flux $Q_{{\mathcal{B}}_{k}}$ is also a prescribed input parameter, and the centroid of an object can be easily computed (see e.g. Zhang & Chen Reference Zhang and Chen2001), we have therefore obtained an expression of the BI equation that is actually usable in practice.

2.2.5 The full boundary integral equation and some remarks

The BI equation (2.10) thus becomes

with the novel bubble term in the last line. The single- as well as the double-layer integrals must be evaluated with the appropriate Green’s functions from (2.14) and (2.18), respectively.

We first remark that (2.33) reduces to the case of the infinite system for $V_{\unicode[STIX]{x1D6E4}}\rightarrow \infty$ (i.e.  $1/V_{\unicode[STIX]{x1D6E4}}=0$ ), as the flux terms vanish and the Green’s functions converge to the well-known expressions for an infinite system, i.e. (Pozrikidis Reference Pozrikidis1992)

for the Stokeslet, and

for the stresslet, where $\hat{\boldsymbol{x}}:=\boldsymbol{x}-\boldsymbol{x}_{0}$ .

Second, the imposed average flow $\langle u_{j}\rangle _{\unicode[STIX]{x1D6E4}}$ can be interpreted as the flow that would prevail in the absence of any objects, and is the most convenient quantity to prescribe a certain flow. A relationship to the corresponding pressure gradient is easily established (Zhao et al. Reference Zhao, Isfahani, Olson and Freund2010, equation (8)).

Third, the prescribed fluxes for the bubbles cannot be chosen arbitrarily. To see this, consider on the one hand

where we used once again the periodicity of the velocity $\boldsymbol{u}$ and the opposite signs of the normal vector $\boldsymbol{n}$ on opposing sides of the unit cell surface $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ . On the other hand,

where the divergence theorem and the incompressibility of the velocity in $\unicode[STIX]{x1D6FA}$ and ${\mathcal{C}}$ together with the definition of the flux have been employed. Combining these two equations, we find

Hence, the fluxes must be chosen such that the total flux is zero, i.e. that the outer fluid volume is conserved. This implies that at least two bubbles are required that oscillate out of phase for the periodic system. Furthermore, if the ambient fluid domain $\unicode[STIX]{x1D6FA}$ is not simply connected (imagine a tube separating the unit cell $\unicode[STIX]{x1D6E4}$ into an inner and an outer domain), the fluxes must be chosen such that the volume within the respective domains is conserved.

Fourth, the initial assumption that all objects are completely located within the unit cell can be relaxed. Surface integrals can be evaluated continuously over the objects’ surfaces even if these surfaces cross the boundary of the unit cell, as is common in simulations of dense suspensions. This property is usually silently assumed in the literature, although it is a priori unclear if it holds due to the linear part in the stresslet (2.18) and the non-periodic centroid term in (2.33). With volume-changing objects it actually follows from the non-trivial interplay between the integrals over the unit cell’s surface $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ and the centroid term. We prove this result explicitly in appendix B.

2.3.1 Basic equation

We now assume that all objects have surfaces of Lyapunov type (i.e. they are ‘smooth’: see Kress (Reference Kress2014), p. 96, for more details). If corners or edges within the surfaces existed, the results would change: see e.g. Pozrikidis (Reference Pozrikidis1992, p. 37) and Kress (Reference Kress2014, chap. 2.5). For smooth objects the single-layer potential is continuous (Kohr & Pop Reference Kohr and Pop2004, chap. 3.4.4) if $\boldsymbol{x}_{0}$ is moved across the surface and the double-layer potential makes a jump (Pozrikidis Reference Pozrikidis1992, equation (2.3.12)). Following these two references, we obtain the Fredholm boundary integral (FBI) equation as

Note that the evaluation point $\boldsymbol{x}_{0}$ is located directly on the surfaces of the objects. The single- as well as the double-layer integrals exist as improper integrals (Kress Reference Kress2014, chap. 6). Equation (2.38) corresponds to the first two lines in (2.1). The missing two terms will be introduced in § 2.3.3 to ensure uniqueness of the solution.

The above FBI equation can in principle be used to find the unknown quantity – if the solution were unique in all cases. In our application presented in § 3, we prescribe the traction/traction-jump $\boldsymbol{F}$ on all objects (also on walls for efficiency reasons). This leads to a so-called Fredholm equation of the second kind that is amenable to the Fredholm theory. As will be analysed and fixed below, equation (2.38) has multiple solutions if bubbles are included. Without bubbles, the solution is unique, as will also be shown below.

On the other hand, prescribing the velocities on all objects leads to a Fredholm equation of the first kind which has various unfavourable properties: the solution is in general not unique, the condition number grows with resolution (compare Karrila & Kim Reference Karrila and Kim1989, p. 127, and Marin, Gustavsson & Tornberg Reference Marin, Gustavsson and Tornberg2012) and no general mathematical theory exists. These are the reasons why alternative approaches for such problem statements have been invented, e.g. the completed double-layer boundary integral method (Karrila & Kim Reference Karrila and Kim1989; Kohr & Pop Reference Kohr and Pop2004; Zhao et al. Reference Zhao, Shaqfeh and Narsimhan2012).

Finally, prescribing the velocities on some objects and the tractions on others yields a mixed type equation. Similar to the first type, no general theory exists and at least parts of the system have ‘difficult’ properties.

2.3.2 Fredholm theory and the non-uniqueness of solution

Henceforth, we consider the case when (2.38) is a Fredholm integral equation of the second kind, i.e. when all velocities are unknown. In order to apply the Fredholm theory, we need to introduce the homogeneous version of (2.38),

where $\boldsymbol{h}$ denotes an eigensolution to the eigenvalue $1$ . Note again that the double-layer integral is meant to be absent for  objects, especially walls. The corresponding adjoint equation (Pozrikidis Reference Pozrikidis1992, p. 106, and Kress Reference Kress2014) is given by

with the eigensolution $\boldsymbol{a}$ and the abbreviation

The integral kernels and their adjoints are weakly singular (see Kress Reference Kress2014, p. 31 and theorem 4.12, Pozrikidis Reference Pozrikidis1992, pp. 36 and 113, as well as Karrila & Kim Reference Karrila and Kim1989, p. 137). This means that all occurring integral operators are compact (Kress Reference Kress2014, theorem 2.30), and that the eigensolutions of the homogeneous and adjoint equations are continuous (see Kress Reference Kress2014, p. 58).

For walls (i.e. open objects) we adopt the convention that closed surface integrals  $\oint$ go over both sides. Due to the continuity of the eigensolutions, however, they provide no contribution. Alternatively, as in the derivation of the BI equation, one can also revert back to walls with finite thickness and take the limit afterwards. The formulas in the following have to be interpreted in the same way. Notice that in the adjoint equation integrals over objects appear which actually have . The compactness of the integral operators also implies that the Fredholm theory can be used to make precise statements about uniqueness and existence of solutions (see Pozrikidis Reference Pozrikidis1992, p. 114, and Kress Reference Kress2014, p. 55f.). For the present purpose the major theorem can be written as follows.

To arrive at uniqueness and existence statements therefore requires us to know all solutions of the adjoint equation. In the case of (2.47), the solutions are

That these are indeed solutions can be easily shown with the help of (2.21). To show that they are the only solutions requires a somewhat longer procedure, similar to § 2.3.4 (we skip it as it is not of any major interest here). Thus, the homogeneous equation (2.39) also has $N_{{\mathcal{B}}}$ solutions. Furthermore, we have, and all solutions $\boldsymbol{a}^{(m)}$ satisfy condition (2.42) due to equation (2.22). Hence, by virtue of the Fredholm alternative, the FBI equation (2.38) has more than one solution if bubbles are included.

2.3.3 Ensuring uniqueness: the full equation

Equation (2.38) does not have a unique solution because the flux of the bubbles is not determined by the equation. To remove this ambiguity, we introduce additional terms into equation (2.38) in such a way that the solution of the new equation is unique and simultaneously also a solution of the old equation (2.38). In analogy to Nie et al. (Reference Nie, Tanveer, Dupont, Li, Xiaobing and Schulze2002), the modified equation is then given by (2.1), which is repeated here for convenience:

This is the central equation that is solved in our VCO-BIM. It is a direct method as the solution $\boldsymbol{u}$ is the physical velocity rather than an auxiliary field. The fact that (2.44) has exactly one solution will be proved below and constitutes a central result of the present work; $\boldsymbol{z}^{(k)}$ can be chosen arbitrarily, as long as the restriction

is satisfied. A convenient choice is

where $S_{O_{k}}$ is the surface area of the $k$ th object. Furthermore, $Q_{O_{k}}$ is the prescribed flux of object $O_{k}$ , which must be zero for all entities except the bubbles.

The integral term in the last line of (2.44) can be interpreted as part of the so-called Wielandt deflation procedure (Kim & Karrila Reference Kim and Karrila1991) for objects with viscosity ratios . This method is sometimes used to accelerate the convergence rate of iterative solution algorithms (Zinchenko, Rother & Davis Reference Zinchenko, Rother and Davis1997; Zinchenko & Davis Reference Zinchenko and Davis2000), but is otherwise optional for . Choosing not to use it amounts to setting $\boldsymbol{z}^{(k)}=0$ (in which case condition (2.45) must be disregarded). For bubbles  that oscillate ( $Q_{O_{k}}\neq 0$ ), however, an additional new term including the surface flux $Q_{O_{k}}$ needs to be taken into account. Note that the last line is an essential ingredient to ensure uniqueness of the solution for bubbles (with and without volume changes), contrary to the usual situation found in the literature. We also remark that the FBI equation remains valid in an infinite system ( $\unicode[STIX]{x1D6E4}\rightarrow \mathbb{R}^{3}$ ) similar to the BI equation from § 2.2.5.

The solution of the patched equation (2.44) is still a solution of the old equation (2.38). This can be shown by multiplying (2.44) with the normal vector $n_{j}$ , summing over $j$ and integrating over some object’s surface $\unicode[STIX]{x2202}O_{k}$ . Using relations (2.21), (2.22) and (2.45) gives $\oint _{\unicode[STIX]{x2202}O_{k}}u_{i}n_{i}\,\text{d}S=Q_{O_{k}}$ , meaning that the flux out of or into the object matches with the prescribed value $Q_{O_{k}}$ , as expected. Moreover, substituting it back into (2.44) recovers the original equation (2.38).

Despite the patch, equation (2.44) is still a Fredholm integral equation of the second kind for the velocities on all objects if $\boldsymbol{F}$ is given. In order to apply the Fredholm theory, we need to introduce again the adjoint of the homogeneous equation, namely

with the abbreviation $M_{j}$ from (2.41); $\boldsymbol{a}$ denotes again the eigensolutions to the eigenvalue $1$ .

The goal now is to prove that (2.47) has only the obvious solution $\boldsymbol{a}=0$ . Theorem 1 then implies that the actual FBI equation (2.44) has exactly one solution. Unfortunately, the procedure used by Pozrikidis (Reference Pozrikidis1992, p. 116f.) cannot be adapted for the proof in the periodic system because the artificial flow that he introduces is a source field. This works in infinite domains where the fluid can escape to infinity, but violates the conservation of the outer fluid volume in periodic domains (even if all objects were volume conserving). We keep the proof rather general, as none has been published before for the periodic system to the best of our knowledge.

2.3.4 Proof of existence and uniqueness of the solution

As a start, we assume that there is at least one non-trivial solution $\boldsymbol{a}$ of the adjoint equation (2.47). We then define an artificial velocity field similar to Odqvist (Reference Odqvist1930, §4) by

This field has a few special properties. First of all, $\boldsymbol{A}$ is defined within the whole space $\mathbb{R}^{3}$ because it inherits the domain from the periodic Stokeslet and because such a single-layer integral exists in the improper sense if $\boldsymbol{x}_{0}$ is located on any surface. Moreover, because the eigensolutions $\boldsymbol{a}$ of the adjoint equation are continuous as explained above, $\boldsymbol{A}$ is continuous as $\boldsymbol{x}_{0}$ crosses any object surface $\unicode[STIX]{x2202}O$ (Kohr & Pop Reference Kohr and Pop2004, chap. 3.4.4). The field is also periodic due to (2.20a ), and we have

due to equation (2.16) and

because of (2.17).

Furthermore, if we define the associated pressure as

with the Green’s function $\boldsymbol{p}$ for the pressure from (2.15), one can show with the help of the singular Stokes equation (Pozrikidis Reference Pozrikidis1996, equation (2.2)) as well as (2.14) and (2.15) that the Stokes equation with $\boldsymbol{A}$ is satisfied everywhere but on the surfaces, i.e. 

$\{O^{(\unicode[STIX]{x1D736})}\}$ denotes the objects and all of their periodic images. The Stokes equation can alternatively be written as

where the stress tensor is given by $\unicode[STIX]{x1D70E}_{ij}^{\boldsymbol{A}}:=-P^{\boldsymbol{A}}\unicode[STIX]{x1D6FF}_{ij}+\unicode[STIX]{x1D707}(\unicode[STIX]{x2202}A_{i}/\unicode[STIX]{x2202}x_{j}+\unicode[STIX]{x2202}A_{j}/\unicode[STIX]{x2202}x_{i})$ . Continuing, the traction (cf. (2.7)) associated with $\boldsymbol{A}$ at the outside (  $\boldsymbol{f}^{\boldsymbol{A},+}$ ) and inside surface (  $\boldsymbol{f}^{\boldsymbol{A},-}$ , normal vector to the outside) of some object $O_{k}$ can be expressed as (compare Kohr & Pop Reference Kohr and Pop2004, equation (3.4.61), and Odqvist Reference Odqvist1930, equation (2.15))

where $M_{j}$ was defined in (2.41). Solving equations (2.54) for $\boldsymbol{a}$ and $\boldsymbol{M}$ leads to

The last required property of the artificial field $\boldsymbol{A}$ is the energy conservation. Following Pozrikidis (Reference Pozrikidis1992, chap. 1.5) and using (2.49) and (2.52), we can derive

The strain rate tensor is defined as

The integral over the unit cell’s surface $\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}$ is simply zero. This follows similar to the derivation of the double-layer integral in § 2.2.4 because $\boldsymbol{A}$ is periodic, $\boldsymbol{f}^{\boldsymbol{A},+}$ contains at most a linear component (due to the definition of the traction, equation (2.7), the pressure, equation (2.51) and the linear term in the pressure Green’s function, equation (2.15)), $\boldsymbol{A}$ is incompressible according to (2.49) and because the average flow is zero as given by (2.50). Furthermore, a similar equation can be derived for the inside of the objects since $\boldsymbol{A}$ is defined everywhere. In the end, we obtain

and

The inequalities follow because the viscosity $\unicode[STIX]{x1D707}$ is ${>}0$ and the integrals contain kernels that are obviously greater or equal to zero.

With all required properties of the artificial field $\boldsymbol{A}$ established, we now proceed to show that the adjoint equation (2.47) does not have any non-trivial solution $\boldsymbol{a}$ . We will do this by a reductio ad absurdum argument. Hence, assume that there is at least one non-trivial solution denoted by $\boldsymbol{a}$ . Following Odqvist (Reference Odqvist1930, § 4), we begin by substituting (2.55) into the adjoint (2.47), leading to

We now multiply with $A_{j}$ , sum over $j$ and integrate over the surface of all objects. The contribution from the second term is simply zero because we can use (2.48) and write

Thus we find

The inequality signs follow from the energy conservation (2.58) and (2.13). Both inequality signs together imply

and due to (2.58a ) we thus have $E_{ij}^{\boldsymbol{A}}(\boldsymbol{x})=0$ for all $\boldsymbol{x}\in \unicode[STIX]{x1D6FA}$ and $i,j=1,2,3$ . This in turn means that $\boldsymbol{A}$ can only represent rigid-body motion within $\unicode[STIX]{x1D6FA}$ (Pozrikidis Reference Pozrikidis1992, chap. 1.5), i.e. 

$\boldsymbol{U}^{\unicode[STIX]{x1D6FA}}$ and $\unicode[STIX]{x1D74E}^{\unicode[STIX]{x1D6FA}}$ are constants that do not depend on $\boldsymbol{x}$ . The symbol ‘ $\times$ ’ denotes the cross-product. Furthermore, using the periodicity of $\boldsymbol{A}$ , we immediately find $\unicode[STIX]{x1D74E}^{\unicode[STIX]{x1D6FA}}=0$ .

Next, we recall that $\boldsymbol{A}$ is continuous across the objects’ surfaces, i.e.  $\boldsymbol{A}|_{\unicode[STIX]{x2202}O}=\boldsymbol{A}|_{\unicode[STIX]{x1D6FA}}=\boldsymbol{U}^{\unicode[STIX]{x1D6FA}}$ , to derive

Here we used the definition of the traction from (2.7), the divergence theorem, the symmetry of the stress tensor and finally the Stokes equation (2.53). Summing expression (2.58b ) over all objects, using (2.64) and that the integral arguments on the right-hand side of (2.58b ) are positive, we find $E_{ij}^{\boldsymbol{A}}(\boldsymbol{x})=0$ for all $\boldsymbol{x}\in O_{k}$ , $k=1,\ldots ,N_{O}$ and $i,j=1,2,3$ . Hence, $\boldsymbol{A}$ must also represent rigid-body motion within every object:

The $2N_{O}$ constants $\boldsymbol{U}^{(k)}$ and $\unicode[STIX]{x1D74E}^{(k)}$ could in principle be different for each $k$ . But because of the continuity of $\boldsymbol{A}$ across $\unicode[STIX]{x2202}O_{k}$ we have $\boldsymbol{U}^{(k)}=\boldsymbol{U}^{\unicode[STIX]{x1D6FA}}$ and $\unicode[STIX]{x1D74E}^{(k)}=0$ for all $k=1,\ldots ,N_{O}$ . All in all, we derived the explicit expression

for the artificial field $\boldsymbol{A}$ .

Next, we exploit the Stokes equation (2.52) which immediately leads to $\unicode[STIX]{x1D735}P^{\boldsymbol{A}}(\boldsymbol{x})=0$ for $\boldsymbol{x}\in \mathbb{R}^{3}\setminus \unicode[STIX]{x2202}\{O^{(\unicode[STIX]{x1D736})}\}$ . The pressure associated with $\boldsymbol{A}$ is thus a simple constant in every connected set, which we will write as

Only the values within the unit cell will be needed (the pressures within the periodic images might be different at first). The stress tensor is thus

Taking the limit onto the surfaces from the outside and inside and multiplying them with the outer normal vector gives the tractions $\boldsymbol{f}^{\boldsymbol{A},+}$ and $\boldsymbol{f}^{\boldsymbol{A},-}$ . By substituting them into (2.55a ) we obtain

with the constants  for $k=1,\ldots ,N_{O}$ . This result is somewhat remarkable: every solution to the adjoint equation (2.47) must be of the form given by (2.69). It also means that the global linear dependence of the pressure that appears in (2.51) via the Green’s function drops out, which is consistent with expression (2.67).

Now, the initial assumption was that there is a non-trivial solution to the adjoint equation. Since all solutions are of the above form (2.69), there must exist one $k^{\prime }\in \{1,\ldots ,N_{O}\}$ with $\widetilde{C}_{k^{\prime }}\neq 0$ . We thus substitute it into (2.47) for $\boldsymbol{x}_{0}\in \unicode[STIX]{x2202}O_{k^{\prime }}$ , and with the help of (2.21) arrive at

The Wielandt deflation term can be active ( $\boldsymbol{z}^{(k^{\prime })}\neq 0$ ) or inactive ( $\boldsymbol{z}^{(k^{\prime })}=0$ ) for a particular object. If it is active, condition (2.45) and  immediately lead to the contradiction $1=0$ . On the other hand, if the term is inactive, equation (2.70) can only be satisfied for . But this means that by definition $O_{k^{\prime }}$ is a bubble, where we demanded that the Wielandt term is always active. Thus we also get a contradiction. This means that our initial assumption (that there is a non-trivial solution to the adjoint equation) must have been wrong, i.e. (2.47) only has the solution $\boldsymbol{a}\equiv 0$ .

To complete the proof, we use the Fredholm alternative from Theorem 1. The homogeneous equation therefore also has only the trivial solution, and consequently the full FBI equation (2.38) has exactly one solution (existence and uniqueness). This holds as long as the Wielandt term exists for objects $O_{k}$ with . Note that for the Wielandt term does not matter concerning uniqueness of the solution, but may be used to accelerate the numerical convergence as remarked before. We also mention that the above procedure carries over to other systems and Green’s functions such as for an infinite domain without any essential changes.

2.4.1 The traction and the constitutive law

As stated in § 2.2 we prescribe a certain outer traction $\boldsymbol{f}^{+}$ on the surface of the bubbles. This is necessary because the introduction of the traction jump $\unicode[STIX]{x0394}\boldsymbol{f}$ (as is done for capsules) in the BI equation would require the application of the Stokes equation at the inside (compare § 2.2.1 and Pozrikidis Reference Pozrikidis1992, pp. 37 and 143) which is not possible because the inside is a compressible fluid with very low viscosity. Due to this very low viscosity, however, we can neglect the shear stress acting on the inside surface of the bubbles and only the inner pressure $P_{{\mathcal{B}}_{k}}$ will be of relevance. Hence, the outer traction is expressible as (Youngren & Acrivos Reference Youngren and Acrivos1976; Power Reference Power1992)

Note that the minus before the pressure comes from the fact that $\boldsymbol{f}^{+}$ represents the force exerted by the fluid on the membrane, and not vice versa. The traction jump $\unicode[STIX]{x0394}\boldsymbol{f}$ must be determined by some constitutive law for the interface, such as the ordinary Young–Laplace equation

$H$ is the mean curvature, taken to be positive for a sphere. This equation is valid for a spatially constant surface tension $\unicode[STIX]{x1D6FE}_{{\mathcal{B}}_{k}}$ , i.e. for interfaces between two immiscible substances. Additional surfactants can lead to a position dependence and non-zero tangential components (Pozrikidis Reference Pozrikidis2001).

2.4.2 Imposing bubble volume changes

We can now prescribe a certain traction jump and an (in general time-dependent) internal pressure to model an oscillating bubble. The traction can then be computed via (2.71) and substituted into the FBI equation (2.1). This should work fine in principle. However, after the substitution we observe that the $P_{{\mathcal{B}}_{k}}\boldsymbol{n}$ term simply drops out due to (2.22), leaving us unable to enforce a certain pressure and thus any volume changes. This deficiency of the FBI equation is because of the fact that the single-layer potential is incapable of producing any flow with sinks or sources (Pozrikidis Reference Pozrikidis1992, chap. 4.1). This in turn originates from the incompressibility (2.16) of the flow produced by the Stokeslet.

Rather than prescribing a certain internal pressure we therefore prescribe a certain flux $Q_{{\mathcal{B}}_{k}}$ . This is easily implemented as the flux appears explicitly in the patched FBI equation (2.1). For the purpose of solving (2.1) we then set $\boldsymbol{f}^{+}=\unicode[STIX]{x0394}\boldsymbol{f}$ .

3.1.1 Surface forces

The coating of microbubbles leads to deviations from the simple coupling between bubble volume and external pressure expected from the expansion/compression of an ideal gas (Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005; de Jong et al. Reference de Jong, Emmer, Chin, Bouakaz, Mastik, Lohse and Versluis2007; Frinking et al. Reference Frinking, Gaud, Brochot and Arditi2010; Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010; Doinikov & Bouakaz Reference Doinikov and Bouakaz2011; Sijl et al. Reference Sijl, Overvelde, Dollet, Garbin, de Jong, Lohse and Versluis2011). There are various models available for varying types of surface coatings that produce such nonlinear responses (e.g. Sarkar et al. Reference Sarkar, Shi, Chatterjee and Forsberg2005; Faez et al. Reference Faez, Emmer, Kooiman, Versluis, van der Steen and de Jong2013). One that incorporates the special properties of phospholipid-coated bubbles (i.e. the size-dependent state of the coating) was presented by Marmottant et al. (Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005), providing a suitable description for their behaviour (Frinking et al. Reference Frinking, Gaud, Brochot and Arditi2010; Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010; Doinikov & Bouakaz Reference Doinikov and Bouakaz2011; Sijl et al. Reference Sijl, Overvelde, Dollet, Garbin, de Jong, Lohse and Versluis2011; Faez et al. Reference Faez, Emmer, Kooiman, Versluis, van der Steen and de Jong2013; Kooiman et al. Reference Kooiman, Vos, Versluis and de Jong2014). The major ingredient is the introduction of an effective surface tension that depends nonlinearly on the bubble’s size. Such a size dependence is the most important aspect of the model for the purpose of margination. Hence, a more elaborate surface model including surface viscosity (e.g. Paul et al. Reference Paul, Katiyar, Sarkar, Chatterjee, Shi and Forsberg2010) is not required here. The relation which we employ can be divided into three major regions (Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005; Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010) as illustrated in figure 2(b).

1. (i) In the high compression regime, the area available per lipid molecule is smaller than its extent, leading to pronounced buckling as observed by ultrahigh-speed imaging (Sijl et al. Reference Sijl, Overvelde, Dollet, Garbin, de Jong, Lohse and Versluis2011). This is modelled by an effective surface tension of $\unicode[STIX]{x1D6FE}_{soft}\approx 0$ below an effective radius $R_{soft}$ (Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010). The bubble is highly deformable in this state (Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005; Rychak et al. Reference Rychak, Lindner, Ley and Klibanov2006).

2. (ii) With increasing radius, a very narrow elastic regime occurs, extending up to a maximal radius $R_{stiff}$ .

3. (iii) Above $R_{stiff}$ the shell ruptures, leaving floating rafts of lipids on the surface (Borden et al. Reference Borden, Kruse, Caskey, Zhao, Dayton and Ferrara2005; Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005). This leads to the very high surface tension $\unicode[STIX]{x1D6FE}_{stiff}$ of a direct air–water interface.

Due to the smallness of the second regime (Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010), we set in the following $R_{soft}=R_{stiff}$ such that the effective surface tension in our case can be written as

where $R_{eff}:=\sqrt[3]{3V/(4\unicode[STIX]{x03C0})}$ is the instantaneous effective radius and $V$ the bubble volume. The traction jump then follows from the Young–Laplace equation (2.72).

The equilibrium radius $R_{0}$ , i.e. the radius assumed when no ultrasound is present, can be located at different positions relative to the transition radius $R_{soft}$ depending on the process of bubble preparation. Most importantly, it was shown that they can be created in the buckled state ( $R_{0}\leqslant R_{soft}$ ) (Borden et al. Reference Borden, Kruse, Caskey, Zhao, Dayton and Ferrara2005; Rychak et al. Reference Rychak, Lindner, Ley and Klibanov2006; Frinking et al. Reference Frinking, Gaud, Brochot and Arditi2010) as desired for safe transport by default within the blood vessels (Lindner et al. Reference Lindner, Song, Jayaweera, Sklenar and Kaul2002).

3.1.2 Modelling the effect of an acoustic source

Because we are interested in the margination behaviour, the exact form of the oscillations is expected to be irrelevant. We therefore prescribe the flux of the $i$ th bubble as $Q_{i}(t)=A_{i}\sin (2\unicode[STIX]{x03C0}ft)$ to model the effect of an external acoustic source. Here, $t$ is the time, $A_{i}$ the flux amplitude and $f$ the frequency. This results in a volume oscillation of

for the $i$ th bubble, with $V_{i}^{(0)}$ being the initial volume at $t=0$ . To conserve the total outer fluid volume as required by (2.37), both bubbles are set to oscillate out of phase, i.e.  $A_{1}=-A_{2}$ . Although not being entirely realistic, it is mandated on a very fundamental level by the periodic boundary conditions and the incompressibility of the ambient fluid. We do not expect that this small restriction affects the validity of the presented conclusions since (i) bubble–bubble interactions are strongly shielded by the RBCs and (ii) margination hinges upon the stiffness variations of the individual bubbles and is therefore independent of the phase of the oscillations.

Continuing, we emphasize that we impose only the instantaneous bubble volume and not a spherical shape. Hence, the bubbles are still deformable, a property which is crucial for the hydrodynamic interaction with the RBCs.

The most important quantity in the present context is the ratio of the stiff to soft duration which we denominate as $\unicode[STIX]{x1D6FF}=T_{+}/T_{-}$ . Here, $T_{+}$ is the time spent in the stiff state (i.e.  $R_{eff}>R_{soft}$ ) and $T_{-}$ the time in the soft state ( $R_{eff}\leqslant R_{soft}$ ). Since margination would trivially be expected for $\unicode[STIX]{x1D6FF}\gg 1$ , we concentrate on $0\leqslant \unicode[STIX]{x1D6FF}\leqslant 1$ in the present study, in agreement with experiments (Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005; de Jong et al. Reference de Jong, Emmer, Chin, Bouakaz, Mastik, Lohse and Versluis2007; Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010; Sijl et al. Reference Sijl, Overvelde, Dollet, Garbin, de Jong, Lohse and Versluis2011). We remark that $\unicode[STIX]{x1D6FF}$ does not depend on the frequency $f$ .

3.1.3 Bubble parameters

In our simulations we set the surface tensions in the soft and stiff state to $\unicode[STIX]{x1D6FE}_{soft}=0.5\unicode[STIX]{x1D705}_{S}$ and $\unicode[STIX]{x1D6FE}_{stiff}=10\unicode[STIX]{x1D705}_{S}$ , respectively. $\unicode[STIX]{x1D705}_{S}$ is the shear modulus of the red blood cells (see below). These choices suitably describe the stiffness of the bubbles relative to RBCs regarding margination while at the same time ensuring numerical stability. Realistic values of $\unicode[STIX]{x1D6FE}_{stiff}=7\times 10^{-2}~\text{N}~\text{m}^{-1}$ and $\unicode[STIX]{x1D6FE}_{soft}\approx 0$ (Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010; Sijl et al. Reference Sijl, Overvelde, Dollet, Garbin, de Jong, Lohse and Versluis2011) would lead to a numerically very stiff problem and consequently require extremely small time steps. The supplementary information (henceforth SI) shows that $\unicode[STIX]{x1D6FE}_{soft}=0.1\unicode[STIX]{x1D705}_{S}$ and $\unicode[STIX]{x1D6FE}_{stiff}=25\unicode[STIX]{x1D705}_{S}$ do not change the results significantly. Furthermore, we fix $R_{min}=1.7~\unicode[STIX]{x03BC}\text{m}$ and $R_{0}=R_{soft}=2~\unicode[STIX]{x03BC}\text{m}$ which are typical values for microbubbles (Borden et al. Reference Borden, Kruse, Caskey, Zhao, Dayton and Ferrara2005; Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010; Kooiman et al. Reference Kooiman, Vos, Versluis and de Jong2014) (using $R_{0}=1~\unicode[STIX]{x03BC}\text{m}$ leaves the results qualitatively unchanged: see the SI). Taking $\unicode[STIX]{x1D6FF}$ as the major control parameter, $R_{max}$ and the amplitudes $A_{i}$ are uniquely determined via the prescribed volume oscillation law (3.2). Assuming an ideal gas within the bubbles and an atmospheric pressure of 100 kPa, a value of $\unicode[STIX]{x1D6FF}=1$ then corresponds to an acoustic pressure amplitude of approximately $P_{A}\approx 45~\text{kPa}$ , in agreement with experimentally used values (Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010).

In most current applications, ultrasound pressure amplitudes and frequencies are in the kilo-pascal and mega-hertz range, respectively (Kooiman et al. Reference Kooiman, Vos, Versluis and de Jong2014; Lammertink et al. Reference Lammertink, Bos, Deckers, Storm, Moonen and Escoffre2015). Such values lead to strong primary and secondary radiation forces (Rychak, Klibanov & Hossack Reference Rychak, Klibanov and Hossack2005; Johnson et al. Reference Johnson, Vormohr, Doinikov, Bouakaz, Shields, López and Dayton2016) making the bubbles agglomerate in a small spot on the vessel wall opposite to the ultrasound transducer (Dayton et al. Reference Dayton, Klibanov, Brandenburger and Ferrara1999; Rychak et al. Reference Rychak, Klibanov, Ley and Hossack2007; Kilroy et al. Reference Kilroy, Klibanov, Wamhoff, Bowles and Hossack2014). This strong localization is highly undesirable for drug delivery applications where a uniform bubble distribution over the entire vessel wall is required. In contrast, we will show below that ultrasound-triggered margination is able to reliably achieve an isotropic distribution if the ultrasound parameters are chosen such that radiation forces become subdominant. For $P_{A}\approx 45~\text{kPa}$ we therefore keep the acoustic frequency at $f=1~\text{kHz}$ in the following. The magnitude of the primary radiation force is then typically of the order of $|\boldsymbol{F}_{rad}|\approx 10^{-15}~\text{N}$ , meaning that it plays only a secondary role as shown explicitly in the SI. We consequently neglect it in what follows. In order to exploit UTM also at higher frequencies, one can reduce $P_{A}$ as exemplified in the SI where we consider $f\leqslant 10~\text{kHz}$ for $P_{A}\approx 6~\text{kPa}$ .

3.2.1 Blood flow constituents

We model the blood flow by explicitly resolving the red blood cells and treating the surrounding blood plasma as a Newtonian fluid (Chien et al. Reference Chien, Usami, Taylor, Lundberg and Gregersen1966). For our simulations we use mostly $15$ RBCs that are distributed randomly within the blood vessel if not noted otherwise. Each RBC has an initial large radius of $R_{RBC}=4~\unicode[STIX]{x03BC}\text{m}$ (Evans & Fung Reference Evans and Fung1972; Freund Reference Freund2014). Continuing, the RBC shear elasticity is modelled via Skalak’s constitutive energy (Skalak et al. Reference Skalak, Tozeren, Zarda and Chien1973) with a shear modulus of $\unicode[STIX]{x1D705}_{S}=5\times 10^{-6}~\text{N}~\text{m}^{-1}$ (Yoon et al. Reference Yoon, Kotar, Yoon and Cicuta2008; Freund Reference Freund2014) and the typical discocyte shape as the reference geometry. This model also includes an area dilatation modulus that is set to $\unicode[STIX]{x1D705}_{A}=10\unicode[STIX]{x1D705}_{S}$ . Furthermore, we additionally introduce an extra surface dilatation energy $E_{a}=(\unicode[STIX]{x1D705}_{a}/2)(S-S_{0})^{2}/S_{0}$ (Krüger Reference Krüger2012) with the corresponding area dilatation modulus $\unicode[STIX]{x1D705}_{a}=10\unicode[STIX]{x1D705}_{S}$ , the instantaneous surface area $S$ and the reference surface area $S_{0}$ . This leads to area deviations of typically ${\lesssim}4\,\%$ . Moreover, bending forces are modelled according to the Canham–Helfrich law (Canham Reference Canham1970; Helfrich Reference Helfrich1973; Guckenberger & Gekle Reference Guckenberger and Gekle2017) with a bending modulus of $\unicode[STIX]{x1D705}_{B}=2\times 10^{-19}~\text{N}~\text{m}$ (Park et al. Reference Park, Best, Badizadegan, Dasari, Feld, Kuriabova, Henle, Levine and Popescu2010; Freund Reference Freund2014) and the spontaneous curvature set to zero. For numerical efficiency we employ the usual approximation that inner and outer viscosities are equal (Krüger Reference Krüger2012; Zhao et al. Reference Zhao, Shaqfeh and Narsimhan2012; Freund Reference Freund2014; Kumar et al. Reference Kumar, Henríquez Rivera and Graham2014), i.e. the viscosity ratio is $\unicode[STIX]{x1D706}_{RBC}=1$ . As a result, any double-layer integrals over the RBC surfaces vanish, and Wielandt deflation terms cannot appear for RBCs (compare § 2.3.3). Nevertheless, both are still present for the bubbles.

The periodic vessel has a length of usually $24~\unicode[STIX]{x03BC}\text{m}$ and a radius of $R_{Vessel}=11~\unicode[STIX]{x03BC}\text{m}$ . Together with the $15$ RBCs this results in a haematocrit of approximately 16 %, a typical value encountered in capillaries (Klitzman & Duling Reference Klitzman and Duling1979; House & Lipowsky Reference House and Lipowsky1987). The larger simulation from figure 2(a) as well as higher haematocrit values lead to the same results which are presented below and in the SI. Furthermore, one possibility for the boundary condition of the vessel wall would be to set its velocity to zero. This, however, leads to a mixed-type Fredholm integral equation. As explained in § 2.3.1, no general mathematical theory exists and this type can be rather performance intensive although it might work in practice. We therefore follow Freund (Reference Freund2007) and fix the wall’s nodes $\boldsymbol{x}_{i}$ via springs to their original position $\boldsymbol{x}_{i}^{(0)}$ , leading to a traction jump of $\unicode[STIX]{x0394}\boldsymbol{f}=\unicode[STIX]{x1D705}_{W}(\boldsymbol{x}_{i}-\boldsymbol{x}_{i}^{(0)})$ , where $\unicode[STIX]{x1D705}_{W}=6.25\times 10^{6}~\text{N}~\text{m}^{-3}$ is the spring constant. Increasing $\unicode[STIX]{x1D705}_{W}$ by a factor of $5$ does not change results qualitatively as shown in the SI. Thus we end up with a Fredholm integral equation of the second kind having exactly one solution as proved in § 2.3.4.

3.2.2 Hydrodynamics

We use our VCO-BIM for 3-D periodic domains as presented in § 2 to solve the Stokes equation. The core of this method is (2.1) which we solve for an imposed average flow chosen such that the maximal flow velocity in the middle of the vessel is roughly $u_{max}\approx 4.7~\text{mm}~\text{s}^{-1}$ , if not noted otherwise. This value matches with physiological flow velocities in capillaries and arterioles (Popel & Johnson Reference Popel and Johnson2005).

The Stokes equation is a good approximation if the Reynolds numbers are much smaller than unity. For the translational motion we find for our system $Re_{T}=2R_{RBC}u_{max}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D707}\approx 0.03\ll 1$ , where $\unicode[STIX]{x1D707}=1.2\times 10^{-3}~\text{kg}~(\text{s}~\text{m})^{-1}$ is the dynamic viscosity of blood plasma (Skalak, Ozkaya & Skalak Reference Skalak, Ozkaya and Skalak1989) and $\unicode[STIX]{x1D70C}\approx 10^{3}~\text{kg}~\text{m}^{-3}$ its density. A different Reynolds number can be defined based on the radial oscillations as $Re_{R}=(2R_{0})^{2}\unicode[STIX]{x1D70C}f/\unicode[STIX]{x1D707}$ . For $f\leqslant 10~\text{kHz}$ (as used in the SI) this results in $Re_{R}<0.07\ll 1$ . We thus conclude that the Stokes equation can faithfully capture the considered RBC and bubble interactions.

3.2.3 Numerical procedure

The general methodology of our numerical implementation was already explained in § 2.5. Here we only mention the remaining aspects that are specific to the present application. The triangle count for the blood vessel is 630 for the $24~\unicode[STIX]{x03BC}\text{m}$ long channel. Rivara’s longest-edge bisection algorithm (Rivara Reference Rivara1984) is used to refine high curvature and close contact regions for the dynamic objects. Hence, the number of triangles varies over time with typical averages of approximately 1500 for the bubbles and 780 for the RBCs. See the SI for some illustrations. Artificial overlapping between the objects within the channel is further suppressed by the introduction of a short-range repulsive potential $E_{Rep}(r_{ij})=[b/(r_{ij}-l_{m})]\text{exp}[l_{c}/(r_{ij}-l_{c})]$ (Noguchi & Gompper Reference Noguchi and Gompper2005; McWhirter, Noguchi & Gompper Reference McWhirter, Noguchi and Gompper2009) with $r_{ij}$ denoting the distance between two nodes (vertices), $l_{m}$ being the minimal possible distance and $l_{c}$ the distance where the potential smoothly drops to zero. We choose $l_{m}=0.01R_{RBC}$ and $l_{c}=0.125R_{RBC}$ , the latter being of the order of the typical edge length of the initial bubble meshes.

The traction jump on the RBCs for the elasticity and dilatation contributions is computed by differentiating the energies with respect to the mesh vertices, as explained by Krüger (Reference Krüger2012) and Guckenberger et al. (Reference Guckenberger, Schraml, Chen, Leonetti and Gekle2016, § 4.2). The repulsive potential is handled in the same way. Bending forces (for the RBCs) and the mean curvature (required for the bubbles, cf. (2.72)) are obtained via Method C as given by Guckenberger et al. (Reference Guckenberger, Schraml, Chen, Leonetti and Gekle2016). Despite being less precise than some alternative methods described in that work, it proved to be more stable than the others.

SPME errors for the computation of the Green’s functions are kept below ${\lesssim}0.01\,\%$ . Increasing the precision by one order of magnitude did not change the results (see the SI). The Gaussian quadrature rules for the integrals use $7$ Gauss points per triangle for the bubbles and the blood vessel and $4$ points for the RBCs. We solve the integral equation via GMRES with a residuum of max.  $10^{-4}$ . Furthermore, the time evolution from (2.73) is obtained by the adaptive Bogacki–Shampine method (Bogacki & Shampine Reference Bogacki and Shampine1989) with the relative tolerance fixed to $10^{-5}$ and the absolute tolerance set to $10^{-4}R_{RBC}$ (Press et al. Reference Press, Teukolsky, Vetterling and Flannery2007). We use the volume rescaling approach for bubbles and RBCs and additionally the hyperplane method for RBCs to handle any artificial volume drift as explained in § 2.5. No special mesh control scheme was necessary for the vessel and the RBCs due to the nature of the prescribed forces, but for bubbles we use (2.74) where the iteration stops once the maximal displacement falls below $10^{-4}R_{RBC}$ . Typical simulations times are in the 1–2 weeks regime on a recent 20 core Intel system.

3.3.1 Microbubbles with constant surface tensions

In order to illustrate the general effect of margination, we first consider the case when the microbubbles are prepared once in the soft state and once in the stiff state. Figure 3(a) shows two simulations without any volume oscillations. The case $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{soft}$ corresponds to coated bubbles that are always in the soft state. Thus, they have a deformability comparable to the RBCs and remain in the centre of the blood stream together with the erythrocytes. On the other hand, setting $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{stiff}$ models pure bubbles that are much stiffer than the RBCs. Hence, they quickly marginate isotropically to the vessel wall. Similar observations are made in figure 3(b) for bubbles oscillating with a frequency of $f=1~\text{kHz}$ while keeping the effective surface tension constant. These results demonstrate that the volume oscillations by themselves do not strongly affect particle migration for the presently chosen parameters.

Figure 3.

Margination behaviour of purely soft and purely stiff non-oscillating (a) and oscillating (b) microbubbles. $\unicode[STIX]{x1D6FE}=\text{const.}$ in all cases. We depict the radial positions of the centroids of two bubbles with constant effective surface tensions in blood flow as a function of time. The surface tensions are set to $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{soft}=0.5\unicode[STIX]{x1D705}_{S}$ (red/orange) and $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{stiff}=10\unicode[STIX]{x1D705}_{S}$ (black/grey). For (b), bubbles oscillate with a frequency of $f=1~\text{kHz}$ , leading to a variation of $R_{eff}$ between $1.7~\unicode[STIX]{x03BC}\text{m}$ and $2.075~\unicode[STIX]{x03BC}\text{m}$ . Curves for different $\unicode[STIX]{x1D6FE}$ constitute independent simulations. The red blood cells, shown as light grey lines, illustrate the cell-free layer between $7~\unicode[STIX]{x03BC}\text{m}$ and the wall. The vessel radius is $11~\unicode[STIX]{x03BC}\text{m}$ , the haematocrit is fixed to 16 % and the maximal flow velocity is $u_{max}\approx 4.7~\text{mm}~\text{s}^{-1}$ . The soft bubbles ( $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{soft}$ ) remain in the centre, whereas the stiff bubbles ( $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{stiff}$ ) show margination.

3.3.2 Lipid-coated microbubbles with radius-dependent surface tension

To demonstrate ultrasound-triggered margination, we consider two lipid-coated bubbles whose shells are modelled with an effective surface tension as described by (3.1) and that are prepared in the soft state ( $\unicode[STIX]{x1D6FF}=1$ , $R_{eff}=R_{soft}$ ). Figure 4(a) depicts the bubbles’ radial trajectories from a simulation where initially no ultrasound is applied. The bubbles are preferably located in the RBC rich core, in agreement with figure 3(a) and experimental observations (Lindner et al. Reference Lindner, Song, Jayaweera, Sklenar and Kaul2002). This allows for secure travel through the vascular system. Once the volume oscillations are activated after around 4.0 s to model an ultrasound source, the bubbles oscillate periodically in stiffness due to the lipid coating that is modelled according to (3.1). This is in contrast to figure 3 where the surface tension remains constant. Most importantly, we observe fast migration towards the vessel wall within less than one second (see the SI for a movie). This time frame corresponds to a travelled distance of less than 4 mm, highlighting the rapidity of the effect. The cause of the fast margination is the lipid shell: as discussed in the previous section, the coating leads to a stiffening during the high-pressure state of the ultrasound signal and a corresponding softening during the low-pressure state (Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005; Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010). As figure 4(a) clearly demonstrates, the overall behaviour is dominated by the stiff stage, as will be further analysed below.

Figure 4.

Ultrasound-triggered margination. (a) Radial positions of the centroids of two microbubbles coated with lipids, modelled according to Marmottant et al.’s law (Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005; Overvelde et al. Reference Overvelde, Garbin, Sijl, Dollet, de Jong, Lohse and Versluis2010). Since the ultrasound is off at the beginning, the bubbles are soft and thus remain in the vessel interior (effective radius $R_{eff}=2~\unicode[STIX]{x03BC}\text{m}$ ). When the oscillations are switched on at ${\approx}4~\text{s}$ , ultrasound-triggered margination leads to rapid migration to the vessel wall. Here, $\unicode[STIX]{x1D6FF}=1$ , i.e. the bubbles are stiff for the first half of the ultrasound period and soft during the second with their effective surface tension varying in the range $\unicode[STIX]{x1D6FE}\in [0.5;10]\unicode[STIX]{x1D705}_{S}$ . The effective radii alternate between $1.7~\unicode[STIX]{x03BC}\text{m}$ and $2.23~\unicode[STIX]{x03BC}\text{m}$ in each period. See the SI for a movie. (b) Average radial positions of the oscillating bubbles and the RBCs for several different values of $\unicode[STIX]{x1D6FF}$ . The rightmost point corresponds to the limit $\unicode[STIX]{x1D6FF}\rightarrow \infty$ (i.e. always stiff) from figure 3(b). Error bars are determined as explained in the main text and the SI.

We continue to demonstrate the robustness of ultrasound-triggered margination by considering microbubbles that are very soft in equilibrium. The bubbles then spend a much longer portion of the ultrasound period in the soft than in the stiff state ( $\unicode[STIX]{x1D6FF}<1$ ). Figure 4(b) depicts the results. The error bars are determined by considering first the minimal and maximal centroid position of all RBCs/bubbles as a function of time for $t>1~\text{s}$ or after definite margination, second a subsequent temporal average and third a weighted average over simulations with different starting configurations. This is similar to the procedure by Müller, Fedosov & Gompper (Reference Müller, Fedosov and Gompper2014) and is explained in more details in the SI. Figure 4(b) shows that the bubbles are still preferably located at the outside of the RBC rich core for $\unicode[STIX]{x1D6FF}<1$ , even for ratios as low as $\unicode[STIX]{x1D6FF}\approx 0.2$ . The margination is completely suppressed only at small values such as $\unicode[STIX]{x1D6FF}\approx 0.1$ where the soft time is around $10$ times longer than the stiff time. The results from the SI for $P_{A}=6~\text{kPa}$ show a transition at $\unicode[STIX]{x1D6FF}\approx 0.3$ indicating that the precise location of the transition depends on the details of the system set-up but nevertheless happens for $\unicode[STIX]{x1D6FF}\ll 1$ . Thus we can conclude that reliable margination is observed if the soft time is at most three times larger than the time in the stiff state ( $\unicode[STIX]{x1D6FF}\gtrsim 0.3$ ).

The effect that small values of $\unicode[STIX]{x1D6FF}$ are sufficient to trigger ultrasound-triggered margination can be understood qualitatively. During the soft state, shearing by the flow and collisions with red blood cells cause deformations of the bubbles. Both are comparably slow processes. During the subsequent stiff stage, however, a high surface tension forces the deformed object back to a spherical shape much more quickly. More quantitatively, the typical relaxation time towards the spherical rest shape in the stiff state can be estimated as $\unicode[STIX]{x1D70F}_{stiff}=2R_{stiff}\unicode[STIX]{x1D707}/\unicode[STIX]{x1D6FE}_{stiff}\approx 0.1~\text{ms}$ . On the other hand, the time required by the flow to deform the bubble away from the spherical shape in the soft state can be estimated by assuming a simple Poiseuille profile with the centre flow velocity $u_{max}=4.7~\text{mm}~\text{s}^{-1}$ . This leads to a shear rate of $s\approx 544~\text{s}^{-1}$ if the bubble is positioned one diameter ( $2R_{soft}$ ) away from the wall. Hence, $\unicode[STIX]{x1D70F}_{deform}=1/s\approx 1.8~\text{ms}$ , which is approximately one order of magnitude larger than the relaxation time scale $\unicode[STIX]{x1D70F}_{stiff}$ in the stiff state. This explains why the latter dominates the margination behaviour. We remark that the ratio $\unicode[STIX]{x1D70F}_{stiff}/\unicode[STIX]{x1D70F}_{deform}$ is formally equivalent to the capillary number $Ca=2R_{stiff}s\unicode[STIX]{x1D707}/\unicode[STIX]{x1D6FE}_{stiff}$ for the bubbles in the stiff state, i.e.  $Ca$ determines their behaviour while being stiff. The interpretation of the ratio $\unicode[STIX]{x1D70F}_{stiff}/\unicode[STIX]{x1D70F}_{deform}$ , however, is fundamentally different here, as it makes a comparison between the two different states rather than making a statement about only one state.

Figure 5.

Average asphericities of oscillating lipid-coated microbubbles. (a) The left figure shows the result for a single microbubble in an infinite shear flow with shear rate $s=544~\text{s}^{-1}$ as a function of $\unicode[STIX]{x1D6FF}$ . At $\unicode[STIX]{x1D6FF}\approx 0.05$ the deformation increases sharply while for $\unicode[STIX]{x1D6FF}>0.1$ almost no deformation is seen. The limiting value $\unicode[STIX]{x1D6FF}\rightarrow \infty$ is shown as a red solid line at the bottom. The surface tension varies in $\unicode[STIX]{x1D6FE}\in [0.5;10]\unicode[STIX]{x1D705}_{S}$ and the oscillation frequency is $f=1~\text{kHz}$ . Error bars indicate the minimal and maximal asphericities during one ultrasound period. The two inset snapshots show the bubbles with their maximal deformation for $\unicode[STIX]{x1D6FF}=0.032$ (left inset) and $\unicode[STIX]{x1D6FF}=1$ (right inset). (b) Asphericity of the bubbles from figure 4(a) averaged over consecutive time intervals of 0.5 s.

The above argument can be explicitly confirmed by considering an oscillating microbubble in a simple linear shear flow in an infinite domain simulated with VCO-BIM ( $V_{\unicode[STIX]{x1D6E4}}\rightarrow \infty$ , compare § 2.2.5). As a measure of deformation, we extract the asphericity $b:=\left[(\unicode[STIX]{x1D706}_{1}-\unicode[STIX]{x1D706}_{2})^{2}+(\unicode[STIX]{x1D706}_{2}-\unicode[STIX]{x1D706}_{3})^{2}+(\unicode[STIX]{x1D706}_{3}-\unicode[STIX]{x1D706}_{1})^{2}\right]/(2R_{g}^{4})$ from the shape, where $R_{g}^{2}:=\unicode[STIX]{x1D706}_{1}+\unicode[STIX]{x1D706}_{2}+\unicode[STIX]{x1D706}_{3}$ is the squared radius of gyration and $\unicode[STIX]{x1D706}_{1}$ , $\unicode[STIX]{x1D706}_{2}$ and $\unicode[STIX]{x1D706}_{3}$ are the eigenvalues of the gyration tensor (Fedosov et al. Reference Fedosov, Pan, Caswell, Gompper and Karniadakis2011). For reference, the discocyte equilibrium shape of an RBC leads to $b\approx 0.15$ . Figure 5(a) shows that the bubble remains almost spherical for $\unicode[STIX]{x1D6FF}$ values as low as $0.1$ , meaning that the bubble is stiff only during ${\approx}9\,\%$ of the ultrasound period. Only below a rather sharp threshold at $\unicode[STIX]{x1D6FF}\approx 0.05$ the bubble behaves akin to a truly soft object exhibiting strong deformation. Note that this value matches well with the ratio $\unicode[STIX]{x1D70F}_{stiff}/\unicode[STIX]{x1D70F}_{deform}\approx 0.056$ , and reasonably well with the bubbles’ transition from the inner core to the outside in figure 4(b). Furthermore, the value $b\approx 0.015$ at $\unicode[STIX]{x1D6FF}=1$ approximately agrees with the asphericity observed for the simulation in figure 4(a) after the ultrasound was switched on, as depicted in figure 5(b).

3.3.3 Further investigations

Pure margination is an isotropic effect: there is no preferred initial migration direction nor a preferred position close to the wall. We exemplify this in figure 6 where we show the trajectories of migrated bubbles from several simulations representing different system realizations. Obviously, no specific accumulation point exists. This is in contrast to migration induced by buoyancy or radiation forces at much higher frequencies (Dayton et al. Reference Dayton, Klibanov, Brandenburger and Ferrara1999; Rychak et al. Reference Rychak, Klibanov, Ley and Hossack2007; Kilroy et al. Reference Kilroy, Klibanov, Wamhoff, Bowles and Hossack2014; Johnson et al. Reference Johnson, Vormohr, Doinikov, Bouakaz, Shields, López and Dayton2016). Note that the inclusion of radiation forces for the chosen parameters leaves the qualitative results unchanged, as described in detail in the SI. The major influence on figure 6 comes therefore from the initial conditions and the length of the simulations.

Figure 6.

Polar plot of several bubble trajectories (i.e. as viewed from the outlet). The figure shows the $\unicode[STIX]{x1D6FF}=0.74$ and $\unicode[STIX]{x1D6FF}=1$ simulations used for figure 4(b), representing different system realizations. Trajectories only shown for $t>1~\text{s}$ or after definite margination. Rare short-lived migration events to the inside occur. Each bubble in each simulation is shown in a different colour. The outer grey dashed line depicts the vessel radius.

Figure 7(a) shows that increasing the frequency from $f=1~\text{kHz}$ to 10 kHz leaves the qualitative results for the radial position unchanged when radiation forces are neglected. Most interestingly, however, the asphericity is approximately reduced by half in the case of the faster oscillations (figure 7 b). The reason is that for 10 kHz less time within one period is available to deform the bubbles before the stiff state takes over, as suggested by the above time scale estimates. This strongly indicates that higher frequencies reinforce the effect that small values of $\unicode[STIX]{x1D6FF}$ are sufficient to obtain ultrasound-triggered margination. Even more, this serves as a hint that the effect of UTM, which has been overlooked so far, might have provided a noticeable contribution to the effectiveness of microbubbles for targeted drug delivery observed in recent in vivo and clinical studies that used higher frequencies (Lammertink et al. Reference Lammertink, Bos, Deckers, Storm, Moonen and Escoffre2015; Kotopoulis et al. Reference Kotopoulis, Dimcevski, Cormack, Postema, Gjertsen and Gilja2016).

Figure 7.

Influence of frequency: behaviour of two oscillating lipid-coated microbubbles for $\unicode[STIX]{x1D6FF}=1$ and a haematocrit of 16 %, once for a frequency of $f=1~\text{kHz}$ and once for 10 kHz (two distinct simulations); $u_{max}\approx 4.7~\text{mm}~\text{s}^{-1}$ . (a) Radial positions of the centroids. The red blood cells are shown in light grey. (b) Corresponding microbubble asphericities averaged over consecutive time intervals of 50 ms.

Continuing, we demonstrate that UTM intrinsically hinges on the presence of the red blood cells. If they are removed, the result in figure 8(a) is obtained, showing clearly that oscillating lipid-coated microbubbles move towards the centre of the channel for $\unicode[STIX]{x1D6FF}=1$ . This is in notable contrast to figure 4(a), where rapid margination for the same set of parameters is observed. Hence, neglecting the influence of the red blood cells in in vitro experiments can easily lead to conclusions that no longer hold for the in vivo case. On the other hand, with a finite haematocrit, margination is always observed in the sense that the RBCs are located at the interior and the bubbles form an outer layer (figure 8 b). The position of this outer layer depends on the size of the RBC-rich inner region which grows when more RBCs are present (Müller et al. Reference Müller, Fedosov and Gompper2014), with the average radial position being reminiscent of a pitchfork bifurcation. The reason is that the RBCs migrate to the centre and push the bubbles to the outside because the latter are seen as stiffer on average as already described above.

Figure 8.

Influence of haematocrit: (a) only bubbles ( $\text{haematocrit}=0$ ); (b) average radial position for $\unicode[STIX]{x1D6FF}=1$ as a function of haematocrit. Values and error bars extracted from several simulations as in figure 4(b).

Moreover, we depict in figure 9 the influence of the flow velocity. Margination still occurs in all cases, but higher velocities tend to decrease the radial position of the marginated bubbles. The effect, however, is comparably small. The upper horizontal axis displays a corresponding effective non-dimensional shear rate defined here by

with an effective RBC diameter $D_{RBC}:=\sqrt{A_{RBC}/\unicode[STIX]{x03C0}}$ and the RBC surface $A_{RBC}\approx 137~\unicode[STIX]{x03BC}\text{m}^{2}$ . This definition is similar to the one by Müller et al. (Reference Müller, Fedosov and Gompper2014) except that we use the maximal instead of the average flow velocity. Hence, we find qualitative agreement with their results for spherical rigid particles (compare figure 3a in Müller et al. Reference Müller, Fedosov and Gompper2014): margination is only a little affected by the shear rate if it is high enough.

Figure 9.

Influence of the velocity for $\unicode[STIX]{x1D6FF}=1$ and a haematocrit of 16 %: average radial position as a function of the flow velocity. Values and error bars extracted from several simulations as in figure 4(b).

Finally, we note that the effect of UTM does not change significantly if the effective surface tension in the soft state is decreased to $\unicode[STIX]{x1D6FE}_{soft}=0.1\unicode[STIX]{x1D705}_{S}$ or the stiff tension is increased to $\unicode[STIX]{x1D6FE}_{stiff}=25\unicode[STIX]{x1D705}_{S}$ , if the size of the bubbles is halved to $R_{0}=1~\unicode[STIX]{x03BC}\text{m}$ , if the initial particle distribution is varied or if the simulation box length is doubled as shown in the SI. Moreover, the average of the nodal velocity of the vessel wall (§ 3.2.1) after a short initial startup phase is ${\lesssim}0.08u_{max}$ and the nodes move less than $0.1R_{RBC}$ , indicating that the spring wall suitably replaces a completely stationary wall. Nevertheless, we also show in the SI that the results do not change if the vessel wall is made five times stiffer.