A fast algorithm for simulating multiphase flows through periodic geometries of arbitrary shape

A fast algorithm for simulating multiphase flows through periodic geometries of arbitrary shape

Gary R. Marple Department of Mathematics, University of Michigan, Ann Arbor, MI, 48109, USA. email:gmarple@umich.edu.    Alex Barnett Department of Mathematics, Dartmouth College, Hanover, NH, 03755, USA. email:ahb@math.dartmouth.edu.    Adrianna Gillman Computational and Applied Mathematics, Rice University, Houston, TX, 77005, USA. email:adrianna.gillman@rice.edu.    Shravan Veerapaneni Department of Mathematics, University of Michigan, Ann Arbor, MI, 48109, USA. email:shravan@umich.edu.
September 10, 2019
Abstract

This paper presents a new boundary integral equation (BIE) method for simulating particulate and multiphase flows through periodic channels of arbitrary smooth shape in two dimensions. The authors consider a particular system—multiple vesicles suspended in a periodic channel of arbitrary shape—to describe the numerical method and test its performance. Rather than relying on the periodic Green’s function as classical BIE methods do, the method combines the free-space Green’s function with a small auxiliary basis, and imposes periodicity as an extra linear condition. As a result, we can exploit existing free-space solver libraries, quadratures, and fast algorithms, and handle a large number of vesicles in a geometrically complex channel. Spectral accuracy in space is achieved using the periodic trapezoid rule and product quadratures, while a first-order semi-implicit scheme evolves particles by treating the vesicle-channel interactions explicitly. New constraint-correction formulas are introduced that preserve reduced areas of vesicles, independent of the number of time steps taken. By using two types of fast algorithms, (i) the fast multipole method (FMM) for the computation of the vesicle-vesicle and the vesicle-channel hydrodynamic interaction, and (ii) a fast direct solver for the BIE on the fixed channel geometry, the computational cost is reduced to per time step where is the spatial discretization size. Moreover, the direct solver inverts the wall BIE operator at , stores its compressed representation and applies it at every time step to evolve the vesicle positions, leading to dramatic cost savings compared to classical approaches. Numerical experiments illustrate that a simulation with can be evolved in less than a minute per time step on a laptop.

S

tokes flow, periodic geometry, spectral methods, boundary integral equations, fast direct solvers

1 Introduction

Suspensions of rigid and/or deformable particles in viscous fluids flowing through confined geometries are ubiquitous in natural and engineering systems. Examples include drop, bubble, vesicle, swimmer, and red blood cell (RBC) suspensions. Understanding the spatial distribution of such particles in confined flows is crucial in a wide range of applications including targeted drug delivery [30], enhanced oil recovery [44], and microfluidics for cell sorting and separation [31]. In several of these applications, the long-time behavior of the suspension is sought. For example: What is the optimal size and shape of targeted drug carriers that maximizes their ability to reach the vascular walls escaping from flowing RBCs [30, 15]? What is the optimal design of a microfluidic device that differentially separates circulating tumor cells from blood cells [49]? More generally, one is interested in estimating the rheological properties of a given particulate suspension in an applied flow, electric, or magnetic fields. A common mathematical construct that is employed in such a scenario is the periodicity of flow at the inlet and the outlet. Therefore, the natural computational problem that arises is to solve for the transient dynamics of a particulate flow through a confined periodic channel driven either by pressure difference or other stimuli.

In this work, we consider periodization algorithms for vesicle suspensions in confined flows. Vesicles—often considered as mimics for biological cells, especially RBCs [54]—are comprised of bilipid membranes enclosing a viscous fluid and their diameter is typically less than 10. The membrane mechanics is modeled by the Helfrich energy [26] combined with a local inextensibility constraint. At the length scale of the vesicles, the Reynolds number is extremely small, therefore, the Stokes equations are employed to model the fluid interior and exterior to the vesicles. The suspension dynamics of this system is governed by the nonlinear membrane forces, the vesicle-vesicle and vesicle-channel non-local hydrodynamic interactions, and the applied flow boundary conditions.

Figure 1.1: (a) Snapshot from a simulation of 1,005 vesicles flowing through an arbitrary-shaped periodic channel. We used 64 discretization points per vesicle and 32,000 points each for the top and bottom walls. The vesicle-vesicle and vesicle-channel hydrodynamic interactions are computed via the Stokes FMM, and the new fast direct solver (Sec. 4) is used to solve the channel BIEs. We used the close evaluation scheme of [6] for the vesicle-to-vesicle and vesicle-to-channel interactions (but not for channel-to-vesicle interactions due to its inapplicability). This simulation took 52 seconds per time step on a laptop with a 2.4 GHz dual-core Intel Core i5 processor and 8 GB of RAM. (b) Plot of the velocity magnitude (red indicates high and blue indicates low) corresponding to the disturbance field generated by the vesicles (obtained by subtracting the pressure-driven “empty pipe” flow from the total velocity field).

Pioneered by Youngren and Acrivos [32, 33], BIE methods are widely used for particulate and other interfacial flows [47]. Their advantages over grid- and mesh-based methods are well-known: reduction in dimensionality, ease of achieving high-order accuracy, and availability of highly scalable fast algorithms. The classical approach for incorporating periodic boundary conditions within the BIE framework is to replace the free-space Green’s function with one that satisfies the periodicity condition. This can be expressed as an infinite sum of source images. The double-layer potential defined on an open curve which defines one period of a channel wall with lattice vector , for instance, can be written as

(1.1)

where u is the fluid velocity, is the Stokes free-space double-layer kernel, and is the density function defined on . Classical algorithms, such as the Ewald summation [16, 25, 46] for accelerating the body calculation that arises from discretizing (1.1), use a partition of unity to split the discrete sum into rapidly converging sums for the nearby and distant interactions, handling the latter in the spectral domain. The local interactions are in number by construction, leaving distant interactions which can be evaluated accurately at the targets by combining local interpolations onto a regular grid with the fast Fourier transform (FFT), a technique named particle-mesh Ewald [14]. Recently, an accurate variant called spectral Ewald has been developed for particulate flows [39, 1, 55] in periodic geometries.

Although such FFT-based methods are widely used, owing to their ease of implementation, they suffer from several drawbacks. Firstly, the FFT introduces a complexity, and although the constants are rather small for FFT methods, the scalability of communication costs on multicore architectures is suboptimal (see [19] for a detailed discussion). Secondly, the lack of spatial adaptivity makes them somewhat impractical for problems with multi-scale physics. Finally, the “gridding” required is expensive, becoming even more so in three dimensions. One of the main goals of this paper is to introduce a simple alternative algorithm that is , and, since it exploits existing fast algorithms, overcomes many of these limitations which can be quite restrictive for constrained geometries that have local features (see Fig. 1.1).

Synopsis of the new method. The proposed periodizing integral equation formulation is based on the ideas introduced in Barnett–Greengard [5] for quasi-periodic scattering problems. It uses direct free-space summation for the nearest-neighbor periodic images, whereas the flow field due to the distant images is captured using an auxiliary basis comprised of a small number of “proxy” sources. Periodic boundary conditions are imposed in an extended linear system (ELS) that determines both the wall layer densities (to enforce no-slip boundary condition on the channel) and the proxy source strengths. Although one block of this ELS is rectangular and ill-conditioned, its pseudo-inverse is rapid to compute, allowing accuracy close to machine precision [12, 41]. The disturbance velocities and the hydrodynamic stresses due to the presence of vesicles enter the right-hand side of the ELS so that the combined flow field is periodic from channel inlet to outlet and vanishes on the walls. The proposed integral formulation is versatile in handling the imposed flow boundary conditions: applied pressure-drop across the channel, or imposed slip on the channel (e.g., to model electroosmotic flows), simply modify the right-hand side of the ELS. The scheme can handle various dimensions of periodicity and easily extends to 3D [41].

The proposed particulate flow solver is based on the work of Veerapaneni et al. [53]. It employs a semi-implicit time-stepping scheme to overcome the numerical stiffness associated with the integro-differential equations governing the vesicle evolution. A spectral (Fourier) basis is used to represent the vesicle and channel boundaries. The required spatial derivatives are computed via spectral differentiation and the singular integrals are also computed with spectral accuracy using the product quadrature rule given in Kress [34]. The FMM is used to accelerate the computation of the vesicle-vesicle hydrodynamic interactions. A simple correction term is introduced in the local inextensibility constraint applied at every time step. It eliminates error accumulation over long time-periods that usually leads to numerical instabilities.

Since the channel geometry remains fixed, its BIE linear system may be inverted once and for all. We use a direct solver to precompute its inverse; the channel wall densities can then be determined at each time step for the small cost of a matrix-vector multiply. Since the matrix associated with wall-wall interactions has off-diagonal low-rank structure, the use of a hierarchical fast direct solver reduces the cost of both the precomputation, and of a solve with each new right-hand side to , where is the number of points on the channel walls. Crucially, the cost involved with each new solve is very small when compared against the standard combination of an iterative solver and an FMM: for example, Table 3 shows that, in our setting, the former is three orders of magnitude faster than the latter. The fast direct solver for the ELS is a Stokes version of the periodic Helmholtz scattering solver of Gillman–Barnett [20], with the extra complication that the channel walls form a continuous interface (in [20] it was possible to isolate the obstacles). The continuous interface demands a new compression scheme, presented in Sec. 4. While fast direct solvers have received much attention recently, the authors believe that this work is the first to apply them to particulate flows.

Advantages. The conspicuous advantage of our method is the simplicity that comes from the use of free-space kernels (as opposed to lattice sums or particle-mesh Ewald) and a pressure drop condition that is applied directly. More importantly, this feature allows us to use state-of-the-art high-order quadratures for the singular [36] and nearly-singular integrals [27, 6] and open-source FMM implementations [22]. The periodization scheme itself is spectrally accurate in terms of the number of proxy points. Numerical experiments illustrate that the same number of proxy points are required to achieve a specified accuracy independent of the complexity of the geometry.

Limitations. In this paper, we restrict our attention to two-dimensional problems. Although the periodization scheme itself is straightforward to extend to three dimensions, other components of the numerical method, such as the quadratures and the direct solver for surfaces in three dimensions, require more work. We assume that there are no holes in the given periodic domain and that the viscosities of the fluids interior and exterior to vesicle membranes are the same. Both of these assumptions are merely for simplicity of exposition (e.g., see [48] to relax the latter assumption and the completed double-layer formulation [45] for the former).

Our recent algorithm to evaluate the nearly-singular integrals with spectral accuracy [6] is applied to evaluate the vesicle-vesicle hydrodynamic interactions when they are close to each other. However, we do not apply any close evaluation schemes for the channel-to-vesicle interactions ([6] cannot be applied directly to this setting). This limitation prohibits us from performing simulations of tightly-packed suspensions. One possible remedy is to switch the discretization of the channel from a global basis (Fourier) to local chunk-based schemes, for which close evaluation schemes for Stokes potentials already exist [43].

1.1 Related Work

There is a large body of literature on BIE methods for periodic Stokesian flows of rigid and deformable bodies. A few examples include [60, 38, 17, 23] for fixed particles as in a porous medium, [42, 61, 18, 58, 1] for particulate suspensions, and particularly [59, 51] for vesicle flows. However, almost all studies focus on flows through either simple geometries, such as flat channels or cylinders, or without any constraining walls (i.e., one-, two-, or three-periodic systems in free-space). The work of Greengard–Kropinski [23] uses an intrinsically two-dimensional complex-variable formulation to periodize a FMM-based solver for a fixed doubly-periodic geometry in cost. However, the imposition of pressure-drop conditions in their scheme is a subtle matter involving non-convergent lattice sums. In contrast, in the present work, such conditions are applied simply and directly and the cost remains . The only work that we are aware of where vesicles inhabit an arbitrary periodic geometry is that by Zhao et al [58] where capsules flowing through deformed cylinders were simulated. The constraining geometry is embedded in a box and a Green’s function that satisfies periodic conditions on the box is used. One of the drawbacks of this approach is that a pressure drop cannot be imposed directly but is determined from the mean flow. Furthermore, a large amount of auxiliary data might be introduced because of the embedding in the case of geometries that have multi-scale spatial features. The authors would like to point out that to date, there are no reported results on particulate flows through complex periodic geometries such as those shown in Fig. 1.1.

Fast direct solvers, such as -matrix, HBS, HSS, and HODLR ([57, 50, 56, 8, 9, 21, 2]) which utilized hierarchical low-rank compression of off-diagonal blocks, are naturally applicable to solving the linear systems arising from the discretization of boundary integral equations, thanks to the smoothly decaying property of Green’s functions for far interactions. Since the construction of the fast direct solver dominates the computational cost, it is not beneficial to build a solver for the evolving geometries. Instead, the proposed method constructs a fast direct solver for the fixed constrained walls, then reuses this precomputed solver at each time step to evolve the vesicles. The computational cost for both steps scales linearly with the number of discretization points. The cost for each time step is much reduced compared to an iterative FMM-based solve. While this manuscript describes an HBS solver (see [21] or Section 4) for simplicity of presentation, alternative inversion techniques can be seamlessly substituted in.

1.2 Outline

This manuscript begins by describing the periodic solver for the steady Stokes equation in a channel given velocity boundary conditions (Section 2). The numerical scheme for the evolution of vesicles and the treatment of the coupling between the vesicle and channel BIEs is presented in Section 3. The fast ELS solver for finding the channel densities and the proxy strengths is presented in Section 4. The accuracy, stability, and computational complexity of the method will be illustrated via numerical experiments in Section 5. Finally, the manuscript concludes with a summary and statement of future work in Section 6.

2 Periodization scheme

2.1 Preliminaries: Stokes potentials and the non-periodic BVP

We first define the standard kernels and boundary integral operators used [37] (for our 2D case see [29, Sec. 2.2, 2.3]). Let be the fluid viscosity, a scalar constant. Let be the velocity field and the scalar pressure field for . The pair is a solution to the Stokes equations if

(2.1)
(2.2)

These express force balance and incompressibility, respectively.

The Stokes single-layer kernel (stokeslet) from source point to target point has tensor components

(2.3)

where , , and is the Kronecker delta. Given a density (vector function) on a source curve , the single-layer representation for velocity is then , i.e.,

(2.4)

The associated pressure function is

(2.5)

For the double-layer velocity representation , we have, using the surface normal at each point on the source curve ,

(2.6)

We write its associated pressure function as

(2.7)

We use the notation to indicate the double-layer boundary integral operator from source curve to target , i.e., . If the target and source curves are the same () then is to be taken in the principal value sense and has a smooth kernel for smooth . We have, for a -smooth curve and any , the jump relation

(2.8)

for the interior limit of velocity. Here, is the identity tensor. The non-periodic prototype BVP that we will need to solve is that the pair satisfies the Stokes equations in a bounded domain for given velocity (Dirichlet) data on its boundary . To solve this problem, we insert the double-layer representation into (2.8) to get the 2nd-kind boundary integral equation (BIE) on :

(2.9)

This BVP, and the resulting BIE, has one consistency condition, , and null-space of dimension one corresponding to adding a constant to [29].

Since we will also need to impose traction (Neumann) matching conditions, we need the traction on a target curve due to the above representations. Given a function pair , the Cauchy stress tensor at any point has entries

(2.10)

The hydrodynamic traction of this pair, i.e., the force vector per unit length applied to the fluid at a surface point with outward unit normal , has components

(2.11)

where summation over is implied. Applying (2.11) to the pair, due to the single-layer velocity (2.3) and pressure (2.5) kernel (with fixed source point ), gives the single-layer traction kernel

(2.12)

which we abbreviate by . Likewise, applying (2.11) to the double-layer pair (2.6) and (2.7) gives, after a somewhat involved calculation (eg [40, (5.27)]), the double-layer traction kernel tensor

(2.13)

where we defined the target and source “dipole functions”

respectively. The use of the symbol to mean the traction operator and the double-layer traction kernel will be clear by context. The hypersingular boundary integral operator for the traction of the double-layer from source curve to target we call . To clarify,

where in the final expression, the kernel has tensor components (2.13).

Figure 2.1: Geometry for periodization scheme of Section 2. The grey shows the infinite periodic pipe, and the blue dots the quadrature nodes for the central domain and its near neighbors. The central domain is bounded by , and side walls and , and has the normal senses shown. The proxy points (red) lie on circle of radius .

2.2 Dirichlet problem in a driven periodic pipe

We consider a single unit cell confined by one period of the walls above and below. The full periodic pipe domain is then , where is the lattice vector with period . See Fig. 2.1.

It is conceptually simplest to begin with the following strictly-periodic “empty pipe” problem. Given periodic velocity (Dirichet) data and on the up and down walls, find a solution in that is periodic up to a constant pressure driving per period, i.e.,

(2.14)
(2.15)
(2.16)
(2.17)
(2.18)

The consistency condition on the data is , and the nullity 1, as in the non-periodic case. In our application, the data , will be (minus) the flow velocity induced by a periodized set of vesicles inside .

A standard approach for solving this BVP in the strictly-periodic case would be to sum the double-layer kernel over all periodic copies in order to obtain the periodized version of the kernel:

(2.19)

The representation is then, using to indicate the one period of the upper and lower walls,

(2.20)

By analogy with (2.9), the density could be determined by solving the 2nd-kind integral equation

(2.21)

with .

Remark 2.1

If , this approach could also be used after substracting from velocity data from the Poiseuille flow , , with . The result is a strictly-periodic BVP with and modified data . However, we will find the following approach much more convenient.

Instead, we reformulate the BVP on the single unit cell , introducing a left side wall and right side wall (Fig. 2.1). Note that, given a periodic pipe , the choice of where to place the wall to subdivide the unit cell is arbitrary. We choose them to be vertical for convenience. Furthermore we relax the periodicity condition on the wall velocity data , , and impose between and periodicity conditions for velocity and traction with given arbitrary mismatch , , that we call the “discrepancies” [5]. Thus,

(2.22)
(2.23)
(2.24)
(2.25)
(2.26)

By unique continuation from Cauchy data, if and are periodic, and we choose and , where here indicates the normal on the and walls, then the above BVP is equivalent to (2.14)–(2.18). The special case gives pressure-driven flow in a periodic pipe, free of vesicles. In the general case, there is still a consistency condition on the data. The advantage of the above (non-periodic) BVP in the single unit cell is that the data may be induced by a sum over vesicles which includes only the nearest images and that pressure driving is incorporated naturally.

To solve (2.22)–(2.26), we use a kernel containing only the near-field images, plus a small auxiliary basis for smooth Stokes solutions in to account for the effect of the infinite number of far-field images. For the latter, we use the “method of fundamental solutions” (MFS) basis [7, 4] of stokeslets with sources lying on a circular contour enclosing . (These are also known as “proxy points” [21].) To be precise, the velocity representation is

(2.27)

where

(2.28)

is a sum over free-space kernels living on the walls in the central unit cell and its two near neighbors. The second term contains basis functions that satisfy the Stokes equation in the physical domain living in the central unit cell. The basis needs to accurately represent any field due to the “far” periodic copies (i.e. those indexed ). The source points are equispaced on a circle of sufficiently large radius centered on the central unit cell, and

(2.29)

is the stokeslet at the th proxy source. Each coefficient , totalling unknowns. This may be viewed as approximating a single-layer density lying on the circle which is able to represent inside any field due to sources lying outside. Since the sources are distant from , the convergence is exponential with a rapid rate; we only need (typically less than ) independent of the complexity of the domain or the number of quadrature points needed to accurately represent it.

The pressure representation corresponding to (2.27) is (summing (2.7) in the same fashion),

(2.30)

Thus, for any density and coefficients , solves the Stokes equations in .

Remark 2.2

There are constraints on the radius : larger allows for more rapid error convergence with respect to , but if is larger than , then the circle encloses some image sources and the size of the coefficients grow exponentially large, resulting in catastrophic cancellation. Hence, we fix in this study.

Constructing a linear system is now simply a matter of inserting the representation (2.27) into each of the conditions (2.23)–(2.26), which we now do. Imposing the velocity data on and using the jump relation (2.8) (which only affects the term in (2.28)) gives two coupled boundary integral-algebraic equations,

(2.31)
(2.32)

which we may summarize as

Imposing periodicity in matching velocity and traction data (2.25)–(2.26) gives, after noticing cancellations of all of the close wall-wall interactions,

(2.33)
(2.34)

which we summarize as

Thus, the four coupled boundary integral-algebraic equations (2.31)–(2.34) may be stacked in pairs and written in a block form

(2.35)

The roles of the block matrices are as follows: is a 2nd-kind operator mapping wall densities to (, ) wall velocities, maps auxiliary coefficients to wall velocities, maps wall densities to their discrepancies in the periodicity conditions, and maps auxiliary coefficients to their discrepancies. The functions and contain the Dirichlet and discrepancy data.

2.3 Discretization of the linear system

The coupled integral-algebraic equations (2.35) need to be discretized: this is performed in a standard fashion using quadrature nodes on , nodes on , and nodes on (the nodes on being those on displaced by ). The nodes on and are generated using the periodic trapezoid rule applied to a smooth parametrization of the curves, while the and wall “collocation” nodes are chosen to be Gauss–Legendre in the vertical coordinate (no weights are needed for these nodes) [35, Ch. 9].

The Nyström method [36, Sec. 12.2] is used to discretize . For instance, given the quadrature weights on , the matrix discretization of the tensor block of has elements

(2.36)

which uses the diagonal limit , where is the curvature of the boundary, and the unit tangent vector. Other blocks are filled similarly. The result is to replace (2.35) by a discrete linear system of similar structure, which will be solved with a fast direct solver described in Section 4. For more details on a similar periodic discretization scheme, see [12].

3 Application to particulate flows

Now we describe the application of the above periodic BVP solution to vesicle flow simulations, extending the scheme of Veerapaneni et al [53] to periodic flow of vesicle suspensions in rigid pipe-like geometries. We begin with just a single (periodized) vesicle with its boundary lying in the periodic unit cell . We solve the steady-state flow problem given forces on the vesicle and then use this solver within the time-stepping scheme.

3.1 Solving for quasi-static fluid flow given the interfacial forces

For simplicity, we consider the case without viscosity contrast. Let parametrize the vesicle membrane according to arc-length . The membrane generates forces on the fluid due to bending and tension , where is the bending modulus and is the tension. The total force at each point on is then . Stress balance and no-slip condition at the membrane-fluid interface imply the jump conditions and respectively, where denotes the jump across . In the case of an isolated vesicle in free-space and assuming is known, the representation for the fluid velocity and the pressure satisfies these jump conditions as well as the Stokes equations in the bulk.

Our goal in this section is, given only the forces on a periodized vesicle , , to solve for the fluid flow which results in the periodic confining geometry , with no-slip boundary conditions and given pressure drop across the channel. The basic idea is to write as a sum of the “imposed” flow that the vesicle would generate in an infinite fluid, plus a “response” flow due to the confining geometry . (This is the same concept as the incident and scattered wave in scattering theory [13].) A standard approach in the case might be to use a periodic imposed flow and for the response to solve the periodic BVP (2.14)–(2.18) with velocity data given by the negative of the imposed flow measured on and . The sum of imposed and response flows then would meet our goal. However, this approach has the disadvantage of relying on periodic Greens functions.

Instead, we use the following representation for the physical flow velocity:

(3.1)

The imposed flow (the first term) involves only the vesicle and its immediate neighbor images, as with (2.28). We define the associated imposed pressure similarly: .

The response is then the solution to the single-unit-cell BVP (2.22)–(2.26) with the following data involving traces of the imposed flow on the walls:

(3.2)
(3.3)
(3.4)
(3.5)

It is simple to check that (3.1) then satisfies no-slip velocity data on and , is periodic, and the pressure representation has the required pressure drop (2.18). Note that, as in the block of the previous section, there is cancellation in the discrepancies and so that even when vesicles come close to, or intersect, or , there are no near-field terms. Effectively, the and walls are “invisible” to the vesicles.

To summarize, the algorithm for solving the static periodic pipe flow problem given vesicle forces and the driving has three main steps:

  1. Evaluate the right-hand side data (3.2)–(3.5), ie

    this will be done with the FMM, except for when the vesicle is close to the wall, in which case a recent spectral close evaluation scheme is used [6].

  2. Solve the rectangular linear system (2.35) for the density and coefficient vector ; this is done with the fast direct solver to be described in Section 4.

  3. Evaluate (being the solution to (2.22)–(2.26)) using the representation (2.27); this will again be done via the FMM to get . It is clear that is linear both in and .

This three-step procedure will become one piece of the following evolution scheme.

3.2 Time-stepping scheme

So far we have only described a quasi-static solution for driven by and . To close the system, we enforce no-slip conditions on the vesicle, , where . Substituting (3.1) gives the first equation in the integro-differential system of evolution equations for the vesicle dynamics, namely

(3.6)
(3.7)

where . Unlike other particulate systems (e.g., drops), the interfacial tension is not known a priori and needs to be determined as part of the solution. It serves as a Lagrange multiplier to enforce the local inextensibility constraint—the second equation (3.7) in this system—that the surface divergence of the membrane velocity is zero. This system is driven by , which could vary in time (we take it as constant in our experiments).

The governing equations (3.6)-(3.7) are numerically stiff owing to the presence of high-order spatial derivatives in the bending force. As shown in [53], explicit time-stepping schemes, such as the forward Euler method, suffer from a third-order constraint on the time step size, rendering them prohibitively expensive for simulating vesicle suspensions. Therefore, we use the semi-implicit scheme formulated in [53] with a few modifications to improve the overall numerical accuracy and stability. Given a time step size and the membrane position and tension at the th time step, , , we evolve to by using a first-order semi-implicit time-stepping scheme on (3.6)-(3.7). For implementational convenience, however, we treat the discretized membrane velocity, denoted with a slight abuse of notation by , as the unknown instead of . The scheme, then, is given by

(3.8)
(3.9)

Since the bending force is a nonlinear function of the membrane position, the standard principle of semi-implicit schemes—to treat the terms with highest-order spatial derivatives implicitly [3]—has been applied to the particular linearization. The tension is treated implicitly and the vesicle-channel interaction, explicitly111When the vesicle is located very close to the channel, say away where is the lowest distance between spatial grid points on the vesicle, a semi-implicit treatment of would be more efficient since the interaction force also induces numerical stiffness in this scenario. Such a scheme, however, requires non-trivial modifications to our fast direct solver of Section 4. Therefore, we postponed this exercise to future work.. The single-layer operator as well as the differential operator are constructed using . In summary, we solve the following linear system for the unknowns :

(3.10)

with given and then update the membrane positions as . The operators are discretized using a spectrally-accurate Nyström method (with periodic Kress corrections for the log singularity [36, Sec. 12.3]) for the single-layer operator and a standard periodic spectral scheme for the differentiation operators. The resulting discrete linear system is solved via GMRES, with all distant interactions applied using the Stokes FMM (e.g., see Appendix D of [53]).

On the initial time step, we generally set . One could obtain an improved initial guess for by setting and solving (3.8) and (3.9) for with . This could be thought of as the first iteration of a Picard iteration for . To enable long-time simulations, we incorporate three supplementary steps in our time-stepping scheme. First, we modify the constraint equation (3.9) as222The main idea here is similar in spirit to the correction formula applied in [52], but in our scheme, we do not introduce any penalty parameters and also rigorously prove its convergence rate (Appendix A).

(3.11)

where represents the perimeter of the vesicle at the th time step. We will refer to this as the “arc length correction (ALC).” We prove in Appendix A that without the ALC, the perimeter of a vesicle will increase monotonically with the number of time steps (total error still scales as ). This would mean that the vesicle’s reduced area can become very low when a large number of time steps are taken. Consequently, the simulated dynamics may correspond to a totally different system than what was originally intended (e.g., a tank-treading vesicle in shear flow might tumble if the reduced area is lowered enough). Executing the ALC at every time step, on the other hand, guarantees a convergence rate, but more importantly, the error is independent of the number of time steps for a fixed (see Theorem A.2). Second, we correct the error incurred in the enclosed area of the vesicle after every time step by solving a quadratic equation in one variable (see Appendix B). Finally, we reparameterize at every time step so that spatial discretization points are located approximately equal arc lengths apart (see Appendix C).

Suspension flow. Although we have presented the case with a single (periodized) vesicle , the above scheme carries over naturally to multiple vesicles. Since such extension has been described previously in other contexts (e.g., see [53] for free-space and [48] for constrained geometry problems) and does not modify in any way our periodization scheme, we only highlight the main steps here. First, the single-layer potential in the representation (3.1) is replaced with a sum of such potentials over all of the vesicles. The equations for the imposed flow data on the walls (3.2)–(3.5) are then modified accordingly. In discretizing the evolution equation for each individual vesicle, the bending force in the self-interaction term is treated semi-implicitly, similar to (3.8)–(3.9). However, the bending forces in the vesicle-vesicle interaction terms can either be treated semi-implicitly, resulting in a dense linear system, or explicitly, resulting in a block tri-diagonal system. While the latter scheme has a marginally smaller computational cost, the former scheme has better stability properties in general since vesicle-vesicle interactions can induce stiffness into the evolution equations when they are located close to each other. In our implementation, we treat all vesicle interactions semi-implicitly. We use a recent single-layer close evaluation scheme [6] to compute the nearby vesicle interactions, whereas for distant interactions, we use the FMM.

4 Fast direct solver for the fixed channel geometry

At every time step, the right-hand side in (3.10) must be evaluated, which involves the channel response 3-step solution given at the end of Section 3.1. However, since the channel geometry is fixed, a fast direct solver enables the second, potentially most expensive, of these three steps to be performed in time with a very small constant. Recall that this step involves solving the block integral equation system (2.35). Upon discretization, this becomes a rectangular system of size given by

(4.1)

This section presents a fast direct solution technique for (4.1). The idea is to precompute for computational cost, the factors in the block matrix solve. Then, the contribution from the channel geometry in the time-stepping scheme only requires a collection of inexpensive linear-scaling matrix vector multiplies.

This section begins by presenting the block solve. The remainder of the section describes how to efficiently build and apply the block solver. The bulk of the novelty in this work lies in the linear scaling technique for representating the interactions of neighboring geometries.

4.1 The block solve

The solution to the rectangular system (4.1) is given by

where , and denotes the pseudoinverse of . is often referred to as a Schur complement matrix.

The Schur complement and need only be formed once, independent of the number of right-hand sides (i.e. time steps). Beyond that, the array need only be computed once per solve.

When a large number of points are needed to discretize the walls due to a complex geometry or small vesicle size, the cost of computing the block solve is dominated by the cost of computing the inverse of . When is large, computing is computationally prohibitive. Fortunately, the matrix has structure which can be exploited to reduce the computational cost of the block solve.

Recall that is the discretization of (2.28) added to . The underlying structure of can exploited first by considering its expanded form given by

(4.2)

where corresponds to the discretization of the self () and neighbor () channel geometry interactions.

Since the majority of the discretization points on the neighboring geometries are well-separated, potential theory states that and are low rank (i.e. has a rank , where ). Thus, each matrix admits a factorization of the form , where and are of size for . Thus,

, where and .

A consequence of utilizing the low rank factorization is that the inverse of can be computed via the Sherman-Morrison-Woodbury formula

(4.3)

Note, only the inverse of and need to be computed. Section 4.3 describes a technique for inverting with a cost that scales linearly with for most wall geometries. Other fast inversion techniques, such as [57, 50, 56, 8, 9, 21, 2], can be utilized in place of the method in Section 4.3. The matrix is in size and is small enough to be inverted rapidly via dense linear algebra.

4.2 Construction of low-rank factorization for neighbor interactions

The cost of constructing the factorizations of and using general linear algebra techniques, such as QR, is and thus would negate the key reduction in asymptotic complexity gained by using fast direct solvers, such as the one described in Section 4.3 to invert . To maintain the optimal asymptotic complexity, we utilize ideas from potential theory, similar to those in [20]. Unlike [20], the neighboring geometries touch the geometry in the unit cell. As a result, a new technique for computing the factorizations is required. For simplicity of presentation, we describe the new low-rank factorization technique for the matrix . The technique is applied in a similar manner to construct the low-rank factorization of . As in [20], an interpolatory decomposition [24, 11] is utilized.

{definition}

The interpolatory decomposition of an matrix that has rank is the factorization

where is a vector of integers with , and is an matrix that contains an identity matrix. Namely, .

First, the upper and lower part of the geometry is partitioned into a collection of sections via dyadic refinement where the rectangular boxes enclosing each segment get smaller as they approach the neighboring cells. Thus, , where is the number of sections is partitioned into. Figure 4.1 illustrates the partitioning of the walls for the channel with reservoirs in Figure 5.1(d) when compressing the interaction with the right neighbor . The refinement is stopped when the smallest box contains no more than a specified number of points . Typically, is a good choice.

Figure 4.1: Illustration of the dyadic partitioning used to compress the interaction of the walls in the unit cell with the right neighbor, i.e. .

For each portion of the boundary that is not touching , consider a circle concentric with a box bounding with radius slightly less than the distance from the center of the box enclosing to the “wall” where and meet. From potential theory, we know that any field generated by sources outside of this circle can be approximated well by placing enough equivalent charges on the circle. In practice, it is possible to place a small number of “proxy” points spaced evenly on the circle. For the experiments in this paper, we found it is sufficient to have proxy points. Figure 4.2 (a) and (b) illustrate the proxy points for and . An interpolatory decomposition is then found for the matrix denoted that characterizes the interactions between and the proxy points. The result is a matrix and index vector . For touching , a collection of proxy points are placed evenly on a circle with a radius 1.75 times the radius of the smallest circle enclosing . All of the points on that lie inside this proxy circle are called near points. Figures 4.2(c) and (d) illustrates the proxy and near points for . An interpolatory decomposition is then formed for the matrix

, where corresponds to the indices of the portion of that are near . The points picked by the interpolatory decomposition are called skeleton points.

Figure 4.2: (a) and (b) illustration of the proxy points ( x) for , . (c) illustrates the proxy points and the near points (thick green line) on for .

We can now assemble the factors by sweeping through the regions . Let . Then the matrix is a block diagonal matrix where each block is a and .

Remark 4.1

The factorization as described does not have optimal rank. Optimal rank can be obtained by recompressing via additional low rank factorizations while assembling and . Depending on the geometry and the computer, it may or may not be beneficial to do the recompression.

4.3 HBS inversion

As stated previously, the dense matrix has structure that we call Hierarchically Block Separable (HBS). The HBS structure allows for an approximation of to be computed rapidly. Loosely speaking, its off-diagonal blocks are low rank. This arises because is the discretization on a curve of an integral operator with smooth kernel (when the walls are not space-filling). This section briefly describes the HBS property and how it can be exploited to rapidly construct an approximate inverse of a matrix. For additional details, see [21]. Note that the HBS property is very similar to the concept of Hierarchically Semi-Separable (HSS) matrices [50, 10].

4.3.1 Block separable

Let be an matrix that is blocked into blocks, each of size .

We say that is “block separable” with “block-rank” if for , there exist matrices and such that each off-diagonal block of admits the factorization

(4.4)

Observe that the columns of must form a basis for the columns of all off-diagonal blocks in row , and analogously, the columns of must form a basis for the rows in all of the off-diagonal blocks in column . When (4.4) holds, the matrix admits a block factorization

(4.5)

where

and

Once the matrix has been put into block separable form, its inverse is given by

(4.6)

where

(4.7)
(4.8)
(4.9)
(4.10)

4.3.2 Hierarchically Block-Separable

Informally speaking, a matrix is Hierarchically Block-Separable (HBS) if it is amenable to a telescoping version of the above block factorization. In other words, in addition to the matrix being block separable, so is once it has been re-blocked to form a matrix with blocks, and one is able to repeat the process in this fashion multiple times.

For example, a “3 level” factorization of is

(4.11)

where the superscript denotes the level.

The HBS representation of an matrix requires to store and to apply to a vector. By recursively applying formula (4.6) to the telescoping factorization, an approximation of the inverse can be computed with computational cost; see [21]. This compressed inverse can be applied to a vector (or a matrix) very rapidly.