An accurate integral equation method for Stokes flow with piecewise smooth boundaries

Two-dimensional Stokes flow through a periodic channel is considered. The channel walls need only be Lipschitz continuous, in other words they are allowed to have corners. Boundary integral methods are an attractive tool for numerically solving the Stokes equations, as the partial differential equation can be reformulated into an integral equation that must be solved only over the boundary of the domain. When the boundary is at least C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^1$$\end{document} smooth, the boundary integral kernel is a compact operator, and traditional Nyström methods can be used to obtain highly accurate solutions. In the case of Lipschitz continuous boundaries, however, obtaining accurate solutions using the standard Nyström method can require high resolution. We adapt a technique known as recursively compressed inverse preconditioning to accurately solve the Stokes equations without requiring any more resolution than is needed to resolve the boundary. Combined with a periodic fast summation method we construct a method that is O(NlogN)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(N\log N)$$\end{document} where N is the number of quadrature points on the boundary. We demonstrate the robustness of this method by extending an existing boundary integral method for viscous drops to handle the movement of drops near corners.


Introduction
The Stokes equations are used to model slowly moving, highly viscous, fluid flow. They can be thought of as the zero Reynolds number limit of the Navier-Stokes equations. Of the many applications of the Stokes equations, they are often used to model particle suspensions including solid particles [6,3], drops [21,26,27], or vesicles [28]. In addition they often describe very well the flow near a solid boundary, and can therefore be used to derive effective slip boundary conditions for problems at higher Reynolds numbers [1,7].
An advantage of using the Stokes equations over the Navier-Stokes equations is that the Stokes equations are linear and elliptic, allowing us to recast them as a boundary integral equation (BIE). BIEs have several nice properties. All the information needed to solve a BIE is confined to the boundary of the domain; this leads to an immediate dimension reduction. The Stokes equations can be represented as a second-kind Fredholm equation [23]. Assuming the boundary is sufficiently smooth, after discretization the condition number of the resulting linear system is independent of the number of discretization points used, meaning that very highly accurate solutions are obtainable. Traditionally one major drawback of BIEs was the need to solve dense linear systems. However, by using efficient iterative solvers such as GMRES [30], combined with fast matrix vector products [9,19], the cost to solve the N × N dense linear system can be reduced to O(N ) or O(N log N ), where N is the number of discretization points on the boundary of the domain.
In this paper, we will be considering wall-bounded, periodic Stokes flow. For particulate flows, such models are useful because they allow for the computation of various time averaged quantities over a relatively small reference cell, without the need to simulate an unfeasibly large domain. BIEs have been successfully applied to such problems [20,27,35].
Solving PDEs on Lipschitz domains to high accuracy everywhere in the domain is in general quite challenging, independent of the numerical method used; see for example the discussion in [8]. When solving boundary integral equations on domains with corners, the standard Nyström method fails to achieve optimal accuracy [5]. While the underlying equation has a unique solution, the layer density defined on the boundary can become weakly singular at corner points, thus reducing the accuracy of regular quadrature rules, such as composite Gauss-Legendre quadrature. One approach to solving this problem is to cluster additional quadrature points near the corners. This of course dramatically increases the number of unknowns, while at the same time the accuracy of such an approach is limited [5,12]. A recent paper [34] demonstrates an approach that automates both the spatial adaptivity and the order of the quadrature (similar to hp adaptivity) to achieve high accuracy for complicated domains containing several corners. Other methods have been developped, some of which involve elegant kernel-dependent custom quadrature rules [31,29]. We will use a kernel-independent method known as recursively compressed inverse preconditioning [12,13]. This method has the advantage of being relatively simple to implement on top of existing code, and does not introduce any additional unknowns to the linear system that must be solved. We will demonstrate the robustness of this method when applied to periodic Stokes flow by solving problems involving moving viscous drops.

Governing Equations
The governing equations are the steady incompressible Stokes equations, where µ is the viscosity of the fluid.
In this paper we will restrict our attention to periodic channel flow as depicted in Figure 1. For boundary conditions, we will prescribe Dirichlet conditions on the velocity, and enforce periodicity in the x 1 direction on the velocity. In addition we impose a constant pressure drop across the reference cell, so that the pressure itself is not periodic, but the gradient of the pressure is, Figure 1: Sketch of a periodic channel with a Lipschitz boundary. The period of the channel is L in the x 1 direction, and the minimum height of the channel is h. The normal vector on the channel walls points into the bulk fluid.

Boundary Integral Formulation
For clarity of exposition, in this section and the next we will present the boundary integral formulation, and the treatment for the corners on a nonperiodic domain. The periodicity will be addressed in Section 5.
Consider the Stokes equations (1) defined inside a Lipschitz domain, along with Dirichlet boundary conditions on the velocity. For the non-periodic case, we will not be concerned with the pressure. The normal vector along the boundary Γ always points into the fluid. See Figure 2 for a sketch of such a domain.
As given in [24] [Chapter 2], a solution to the forced Stokes equations is given by the velocity and pressure pair The stress tensor σ j at a point x is given by the combination of the velocity and pressure, Here and in the remainder of this paper we have used the Einstein summation convention, where summation over repeated indices is implied. The letters i and k are reserved for later use and are therefore not used as indices.
In R 2 , we have that Here we have used r to denote the Euclidean norm of r.
The tensors G j and T j m are called the Stokeslet and stresslet respectively. From the Stokeslet and stresslet we can define the single-and double-layer potentials, S[q](x) and D[q](x), as a convolution of the Stokeslet and stresslet with a density function q(x) defined over Γ. Explicitly, where n is the unit normal vector pointing into the fluid. Following [11], we will express the solution to (1) as a combination of the single-and double-layer potentials where η > 0 is an arbitrary constant which governs the relation between the single and double layer potentials. To obtain a well-conditioned system η cannot be too large; we will always take η = 1.
This equation is valid for x ∈ Ω. To obtain a boundary integral equation (BIE) we will take the limit of (3) as x approaches a point x 0 ∈ Γ. To do this, we will need the limiting values of the single-and double-layer potentials, Applying these limits to (3) and matching it to the boundary condition g yields to BIE If Γ is Lyapunov smooth, then the operator β is a compact operator, with eigenvalues accumulating at zero. In this case (4) can be analyzed using Fredholm theory. In particular the Fredholm alternative applies, and in [11] this is used to demonstrate the existence and uniqueness of solutions. If, as in our case, Γ is only Lipschitz smooth, then β is not compact so Fredholm theory cannot be applied. Furthermore, the double-layer potential involves the normal vector, meaning that the double-layer kernel cannot be evaluated pointwise. Nonetheless, it can be shown that (4) has a unique solution, even when Γ is only Lipschitz continuous [5]. In this case the density function q is singular at the corner points, however (3) can still be evaluated for any x ∈ Ω, and (4) can be evaluated anywhere on Γ except at the corner points.

Numerical Methods
To evaluate (4), we use the Nyström method, as described in [4] [Chapter 4]. Let γ(s), s ∈ [0, 2π] parameterize Γ. The BIE (4) can be written in the abstract form The kernel K in this case is the kernel of β, i,e., We will approximate the integral in (5) using a composite Gauss-Legendre quadrature scheme of n pan panels and n q quadrature points per panel, where N = n pan n q is the total number of quadrature points, and w n is the quadrature weight corresponding to the quadrature point s n . We will then enforce (7) at the quadrature points s m , m = 1, . . . , N to get the linear system Here the point values of the density function q are unknown. When setting up the composite quadrature rule, care must be taken to ensure that each corner on Γ is at the intersection of two panels. In this way we avoid difficulties arising from trying to evaluate the normal vector on a corner, as the Gauss-Legendre quadrature points cluster near but are never located at panel endpoints.

Singular Quadrature
When t = s, the kernel of the double-layer potential T j m (0)n m (s) has a removable singularity. The Stokeslet however must be handled using a specialized quadrature technique. We will use the approach given in [22]. We begin by writing the single-layer potential in complex variables, where we have introduced a slight abuse of notation to denote q(z) as density function q written as a complex number, i. e., z = x 1 + ix 2 and q(z) = q 1 (x) + iq 2 (x). The only term in the complex formulation of the Stokeslet that does not have a finite limit as τ → z is the term with the logarithm. We will consider this integral over a single panel, Γ k . The panel extends from α 1 to α 2 , with α 1 , α 2 ∈ C, and is parameterized by α ∈ [α 1 , α 2 ]. Let α be parameterized as, where t ∈ [−1, 1] and ψ = (α 2 − α 1 )/2. The log term in the integral above can be rewritten in the form, where τ (α z ) = z. Define t z to be such that α(t z ) = α z . Then we have that which allows us to rewrite the log integral in (8) as The first integral can be evaluated using the fact that To evaluate the second integral, we expand the density q(t) as a polynomial series to obtain, allows us to write I in the compact form The integrals in h can all be computed analytically using a recursion relation. The vector c = {c j } nq−1 j=0 can be computed by solving the Vandermonde system where q = {q(t 0 ), . . . , q(t nq−1 )}, and V is the Vandermonde matrix. Thus I can be approximated by where ω(t z ) is the solution to the linear system h(t z ) = V T ω(t z ). The stability of these computations is discussed in [14], here we simply note that the computations are stable. Note that this linear system is independent of q. For now we will need the values of I only at the quadrature points t 0 , . . . , t nq−1 . Since we have rescaled and rotated the panel Γ k to run from -1 to 1, the corrected weights ω(t z ) can be precomputed for each of the quadrature points t 0 , . . . , t nq−1 .
Quadrature points on adjacent panels will also need to be corrected. How many points to correct depends on accuracy considerations, but for n q = 16, numerical results have shown that if all the panels are the same size (in parameter space), then correcting the four closest quadrature points in each adjacent panel is sufficient. Again, these corrections can be precomputed.

Recursively Compressed Inverse Preconditioning
As previously mentioned, in the case of Lipschitz domains, (4) has a unique solution. However, at the corners, the density function q becomes singular. Therefore Gauss-Legendre quadrature fails to accurately integrate the integrands in (4) near the corners. Local panel refinement around the corners is one way to mitigate this issue, however, this approach adds a potentially large number of new unknowns, see Figure 3. In addition, the condition number of the resulting locally refined linear system grows with the number of quadrature points. An alternative approach, recursively compressed inverse preconditioning (RCIP) [12,13], can achieve the same accuracy as local refinement by performing a simple precomputation. As the desired accuracy requirements increase, the precomputation grows linearly with n sub , however both the size and the condition number of the final linear system remain fixed.
The main idea behind RCIP is to transform the density function q in (4) into a transformed densityq that is piecewise smooth everywhere on Γ. Then Gauss-Legendre quadrature can be effectively used. To create this transformation, the operator β is written as where β • is compact away from the corners, and β * describes the corner interactions. We then define the transformed density as where I is the identity operator.
Using the split (9), and the transformed density (10), we can convert (4) to the transformed BIE forq, where R = (I+β * ) −1 . If we assume g is piecewise smooth, it can be immediately seen thatq must also be piecewise smooth. Since β • is a smoothing operator, then β • Rq will be smooth everywhere. Then, by contradiction, in order for (11) to hold, Iq =q must be piecewise smooth. It remains to discretize (11). One possibility for discretization is to use two meshes: a coarse mesh Γ coa on which β is discretized, and a fine mesh Γ fin , on which R is discretized. To get Γ fin from Γ coa we first designate the two panels on either side of corner j to be the subset Γ * j ⊂ Γ coa , and define Γ • to be Γ coa \ ∪ nc j=1 Γ * j . The panels on either side of corner j in Γ * j are then dyadically refined n sub times to get Γ * j,fin . The fine mesh is defined as Γ fin = Γ • ∪ ∪ nc j=1 Γ * j,fin . An example of this is shown in Figure 4. We will define β • and β * to be the operator β restricted to the domains Γ • and Γ * respectively. Note that β • is a compact operator. The operator β • can Figure 4: Panel discretization using two meshes. The boundary Γ is first discretized with a coarse mesh, Γcoa. The two panels on either side of each corner are denoted Γ * , and the remaining panels are denoted Γ • . The panels in Γ * closest to the corners are dyadically refined n sub times to obtain Γ fin . Each panel, regardless of its size, contains nq Gauss-Legendre quadrature points.
be discretized on Γ coa to get the matrix B • . This matrix is equivalent to the discretization of β on Γ coa , but with the entries corresponding to both source and target being outside Γ • set to zero. The operator R could be discretized on Γ fin . This would however not be much use, since it would introduce a large number of unknowns. Instead we will exploit a forward recursion relation to construct R on a sequence of meshes covering larger and larger portions of Γ * j . To define the recursion relation, we will need to provide some definitions of different types of meshes. For each corner j, define a sequence of meshes Γ and Γ * j,n sub = Γ * j . On each Γ * j, we will have a six panel type-b mesh, Γ * j, b , and a four panel type-c mesh Γ * j, c . The type-b and type-c panels will be related as shown in Figure 5. To interpolate between type-b and type-c meshes defined over the same interval, we introduce the prolongation matrix P bc . In addition, defining W b and W c to be a diagonal matrix containing quadrature weights on the type-b and type-c meshes respectively, we can define the weighted prolongation matrix P W bc as . This matrix will be of size 6n q for scalar problems, or 12n q for vector problems in R 2 . Define R j,1 as We then compute the sequence R j,2 , . . . , R j,n sub using the recursion relation where I • b and B • j, b are the identity matrix and B j, b with the entries in the two panels around the corner zeroed out respectively. The operator F{·} zero pads the argument to turn a matrix defined on Γ * j, c into one defined on Γ * j,( +1)b . Once we have the matrices R j,n sub , we can constructR, the discretization of R on Γ coa . Outside of Γ * , β * is zero, so from the definition of R, we obtain thatR must be the identity matrix over Γ • . Over Γ * j , the discretization of R is just R j,n sub .
To handle the log singularity when assembling the matrix B j,1b and B • j, b , the same techniques as described in Section 4.1 can be used. Some bookkeeping is needed to account for the fact that the panels are not equally sized in parameter space. In some cases adjacent panels will be double or half the size, however the modifications for the quadrature weights can still be precomputed for each case.

Fast multiplication
The fully discrete transformed BIE defined on Γ coa is which can be solved using GMRES. To accelerate the required matrix-vector products from O(N 2 ) to something computationally feasible, fast summation methods are a necesity. To do this, we will rewrite B andR as the block matrices It follows that B • can be expressed as Then (13) can be rewritten as The matrix vector productRq can be done directly, sinceR is just the identity matrix with one 8n q × 8n q block for each corner. The remaining matrix vector products BRq and B * Rq can be computed with a fast summation method, reducing the cost from O(N 2 ) to O(N log N ) or O(N ), depending on the choice of fast summation method. We will use an O(N log N ) method that can naturally handle periodic problems, as will be described in Section 5.

Post-processing
After we have computedq, we evaluate the velocity at a point u ∈ Ω by using the regular quadrature rule whereq =Rq, and K j is defined in (6). This quadrature is as accurate as if we had access to the fine density defined on Γ fin . From the definitions ofq and R, it is clear that outside of Γ * ,q and q are equivalent. On Γ * , the transformed densityq can be thought of as a weight corrected density function.
Both the double-layer potential and the single-layer potentials contain integrals that become near-singular when t − s is small, i.e. the numerical errors grow large when evaluating the integrals at points close to any boundary. To handle this, the integrals are treated in the manner described in [27]. In short, the main idea is similar to that explained in Section 4.1 where the density is expanded as a polynomial series, and the integrals computed analytically using recursion relations.

Periodicity
We now address the periodicity in (2). The periodic single-and double-layer potentials, S P [q](x) and D P [q](x) are defined as an infinite sum of the singleand double-layer potentials defined in Section 3, x ∈ Ω, (14a) With these operators we can define the periodic operator β P [q], Note that periodic sum is over p ∈ Z 2 ; in other words we are implying periodicity in both spatial dimensions, even though in our actual problems the periodicity will be only in one direction. The reason for this is that it allows us to exploit a more efficient fast-summation method. To enforce the periodicity in one dimension we will embed the channel, which is only periodic over a length L in the x 1 direction, in a doubly periodic box of size L × L. A similar approach is used in [35] where a periodic flow in one direction is embedded in two dimensional periodic box.
As mentioned in Section 2, a constant pressure drop p 0 −p 1 in the x 1 direction is applied. In order to impose this this we will add on an unknown mean velocity u to the layer formulation (4), Note that · denotes a volume average over a full periodic cell, including the parts outside Ω. This quantity has no physical meaning beyond imposing a pressure gradient. An alternative approach would be to add a predetermined background flow, as done in [27].
The mean pressure gradient ∇p must balance the net effect of the wall friction [35]. From (15) the wall friction force is equal to η Γ q dΓ [11]. The system we have to solve is thus given by The splits for the two dimensional Stokeslet and stresslet, as well as truncation estimates are derived in [26]. Here we list the results.
To compute the infinite sums in (14) we will use the spectral Ewald method [18,19]. In the spectral Ewald method the infinite sums in S P and D P are split into two parts based on a so-called Ewald decomposition: one that decays exponentially fast in real space, and one that decays exponentially fast in Fourier space. The real space sum can be truncated in real space, while the Fourier sum can be truncated in Fourier space. The splits for the two dimensional Stokeslet and stresslet, as well as truncation estimates are derived in [27]. Here we list the results.

Ewald splits
The discretized periodic single-and double layer potentials are given by where the * in the summations over p indicates that we are excluding any p that makes x − x m − Lp = 0.
The discretized single-layer potential can be rewritten as the split The k = 0 mode in the Fourier expansion has been shown to be 0 [25]. The decomposition parameter ξ > 0 is called the Ewald parameter and determines the relative sizes of the real and Fourier parts of u G (x). Note that u G (x) itself is independent of ξ. When applying the Nyström method to solve for the unknown pointwise density values, we do not wish to evaluate singular cases where x −x n +Lp = 0. This contribution to u G (x) can be skipped in the real space sum, but for the Fourier sum we will have to subtract off the limiting value given by where γ is the Euler-Mascheroni constant. For the discretized double-layer potential, we use the split, For the k = 0 mode, the choicê guarantees zero mean-flow though the reference cell [2]. For the stresslet, the limit lim r→0 (T j m (r) − T R j m (r))q n m w n = 0, so no limiting value needs to be subtracted when r = 0.
To compute these sums, it is necessary to truncate them. Fortunately, both these sums decay exponentially fast, and they can be truncated to a desired tolerance following the estimates in [27]. With appropriate scaling of the decomposition parameter, the real space part is computed in O(N ) time and the Fourier space part is accelerated to O(N log N ) using fast Fourier transform. In order to use FFTs, the source points are spread to a uniform grid where the computations for the Fourier space sum are carried out. The result is then gathered from the uniform grid to the target points. The spreading and gathering is done using truncated Gaussians whose shape parameter is selected to minimize the approximation error for a given support. For efficiency, fast Gaussian gridding [10] can be used in both the spreading and the gathering steps. For more details see [27].
A numerical difficulty in the periodic formulation lies in the fact that the Ewald representation of the stresslet is not translation invariant due to the zero mode (16). This means that the submatrices R j are different for each corner, even if the corners have the exact same shape. To avoid roundoff error, for large n sub , we would like to have a local coordinate system for each corner centered at x = 0. To assemble these matrices, we first assemble the translation invariant part of the matrices R j , and then add on (16) only at the end. Note that this is necessary only for quite large n sub .

Numerical Examples
To test the periodization scheme, we can test with an exact solution. Pressure driven pipe flow creates the well known parabolic flow profile in the x 2 direction. The exact solution for flow through a flat channel with top wall at x 2 = 1 and bottom wall at x 2 = 0 is given by Note that this flow is constant and therefore periodic in the x 1 direction. As can be seen in Table 1, using only 4 panels allows us to achieve a relative error of 10 −7 everywhere up to a distance of 10 −3 from the boundary. Another test is to prescribe Dirichlet boundary conditions on the solid walls. We will prescribe a constant shear flow u = x 2 , 0 as boundary conditions on the top and bottom wall. The bottom wall is now only piecewise smooth. To  Table 1: Maximum relative errors in a pipe flow simulation along lines parallel to the top wall. As the target points get closer to the wall, the integrands become more challenging to evaluate accurately. The minimum error we are able to achieve is limited by the tolerances of our GMRES solver and Ewald summation. recover a shear flow in the interior from the BIE solution, we must prescribe no pressure drop from inlet to outlet, i.e., ∇p = 0. We will use this simple example to test the RCIP algorithm for various values of n sub . Figure 6 shows the results. If we use the standard Nyström discretization without RCIP, we get quite large errors everywhere in the domain; the errors are not localized around the corner. For n sub = 20 the RCIP algorithm allows us to achieve 11 digits accuracy everywhere in the domain, except for close to the corners.To gain additional accuracy, it is possible to recover the actual fine density and use it along with the special quadrature described in Section 4.4 [12] [Section 10]. In every case 20 panels on both the top and the bottom wall are used. Without any special treatment, the recovered solution differs quite a bit from the exact solution everywhere. With n sub = 20 the BIE can achieve 11 digits accuracy everywhere except very close to the corners. If additional accuracy is required near the corners, the fine density function can be recovered and used in conjunction with the special quadrature described earlier.
The number of refinements needed to achieve a desired accuracy depends heavily on the geometry. For example Figure 7 shows an example where n sub = 20 achieves only a maximum of 7 digits of accuracy. For wedge shaped corners, the domain segments Γ j, b are self similar for = 1, . . . , n sub and j = 1, . . . , n c . If the operator β is scale invariant, then B • j, b will be independent of , and the recursion relation (12) becomes a fixed point iteration which can be iterated to find R j to a desired tolerance without specifying n sub [12] [Section 12] . Unfortunately, in our case, because the single-layer potential contains a logarithmic term, β is not scale invariant, so we cannot use this idea, and n sub must be specified a priori. A formulation involving just the double-layer potential [23] would not suffer from this drawback and could be formulated as fixed point iteration. We have chosen the formulation in [11] because numerical experiments have shown it to be more stable for simulations of squeezing drops [36].

Adding Viscous Drops
Earlier work [27] has looked at modelling the movement of drops inside confined periodic geometries. We now demonstrate the robustness and usefulness of RCIP by extending the method in [27] to model the movement of drops near sharp interfaces. A drop is a packet of fluid that is does not mix with the fluid surrounding it. Surface tension forces prevent the mixing of the drop with the bulk fluid. Both the fluid inside each drop, and the bulk fluid, satisfy the incompressible Stokes equations, where µ denotes the viscosity inside the region bounded by Γ . For convenience we will define the viscosity ratio λ = µ /µ 0 . As λ increases the drop behaves more like a rigid particle. At the interfaces the velocity is continuous, however in general the surface forces acting on the interface from the inside and the outside of the drop are not equal. We will denote the jump in the normal force across interface as δf (x). This jump is related to the curvature κ (x), and the surface tension σ (x) by [24] [Chapter 5], where e denotes the rate of strain tensor in Ω and ∇ s is the gradient along the interface.
To nondimensionalize this problem, we will define h to be the minimum height of the channel. We would like the maximum velocity of an empty channel to be 1. To do this, we will introduce a pressure scale h ∇p /8, a length scale h, a velocity scale h 2 ∇p /(8µ 0 ), and a surface tension scale σ 0 to obtain the full nondimensional form of the problem, where the Laplace number Lp is the dimensionless quantity given by h 2 ∇p /(8σ 0 ) [17]. The time scale is now h/U = 8µ 0 /(h ∇p ). The Laplace number serves the same purpose as the Capillary number in other drop simulations [27], but for pressure driven flows it is more natural to nondimensionalize according to the pressure gradient as opposed to a maximum velocity. For the remainder of this paper we will assume a constant surface tension, so ∇ s σ = 0. Allowing for non constant surface tension would not impact the RCIP algorithm in any way. In [26,27] chemical surface active agents are allowed to change the surface tension of the drops. A BIE formulation [36] for the problem is given by where S P Γ and D P Γ are the periodic single-and double-layer operators defined on Γ .
Taking the limit as and taking the limit As in Section 5 we close the system by relating the density function q, still defined only on the solid walls Γ 0 , to the (nomdimensionalized) average pressure gradient, After computing u(x) on the drop interfaces, the drops move and deform according to the velocity, To evaluate this system of ODEs we will use the fourth order adaptive time stepper described in [16]. Further details on time stepping, including a way to maintain consistent acrclength spacing as the drop perimeter changes, can be found in [26,27].
It is important to note that since the solid walls are not moving, the RCIP matrices R j , j = 1, . . . , n c need only be computed at the start of the simulation. This means that after this precomputation the actual time needed each time step to solve the linear system is independent of n sub . In fact, as using RCIP improves the accuracy of the simulation, an adaptive time stepper may find it easier to meet a desired tolerance allowing larger time steps to be taken. Therefore the overall computation time may well be lower if RCIP is used.

Investigation of Accuracy
To test our method, we will perform a numerical experiment. Figure 9 shows snapshots of a simulation. We begin by creating a high resolution reference solution with RCIP for λ = 1 and λ = 5. The high reference solution will be discretized using n pan = 140 on both the solid walls and the drop, and n sub = 20. The time stepping tolerance for the high resolution simulation is 10 −10 . t = 0 t = 0.5 t = 1 Figure 9: Snapshots of a drop moving inside a channel containing a corner. A pressure drop is applied from left to right. By varying the number of points on the drop and the wall we can perform a numerical convergence study, the results of which are shown in Figure 10. We have run two simulations, one with λ = 1 (blue drop), and another with λ = 5 (red drop). As expected, the drop with λ = 5 deforms less than the λ = 1 drop. Lp = 1 for all simulations.
To compute an error, we will look at the ∞ difference of a numerical solution with the computed reference solution. As described in [27], when updating the positions of the drops it is advantageous to work on a uniform grid, instead of the panel Gauss-Legendre points. Using a uniform grid also allows us to easily upsample the non-reference solutions using an FFT to obtain a solution at the same discretization points as the reference solution.
As can be seen in Figure 10, without RCIP, as we refine the number of points on both the drop and the wall, the numerical solutions for both λ = 1 and λ = 5 converge very slowly towards the high resolution reference solution. When we use RCIP with n sub = 20, the numerical solutions are much more accurate, around the level of the time stepping tolerance for the low resolution simulations of 10 −8 .  λ = 1, n sub = 0 λ = 1, n sub =20 λ = 5, n sub = 0 λ = 5, n sub =20 Figure 10: Spatial error study with and without RCIP for the simulation shown in Figure 9. Using RCIP gives much lower errors, and these are near the time stepping tolerance of 10 −8 . Lp = 1 for all simulations.

Multiple Drops
To demonstrate the robustness of the code, we will investigate the movement of multiple drops of different sizes confined in a periodic channel containing a narrow constriction. Snapshots of such a simulation are shown in Figure 11. As the snapshots make clear, there is a clear, visible difference in the modelled drops depending on whether or not we use RCIP. Further proof of the necessity of the proper handling of corners is shown in Figure 12. Since the fluid inside the drops is incompressible, the area of each drop should be conserved. Note that this is not enforced explicitly, nor is area conservation used in the criteria for the temporal adaptivity. Without RCIP, the area after each drop has passed one periodic channel length is conserved only up to around 10 −4 . When applying RCIP, it is below 10 −10 , a full factor of 10 6 improvement.   Figure 11 as a function of time. The error eventually plateaus around 10 −4 without RCIP, but remains below 10 −10 when RCIP with n sub = 50 is used.

Conclusion
Domains with sharp corners pose challenges for boundary integral methods. The layer density defined on the boundary becomes weakly singular at the corners, and therefore standard quadrature rules cannot be accurately applied. We have demonstrated that a technique known as Recursively Compressed Inverse Preconditioning (RCIP) can be used to accurately solve the Stokes equations on two-dimensional domains with corners. This method requires a small amount of precomputation, however it does not add any additional unknowns to the problem, nor does it increase the condition number of the linear system. Numerical experiments have demonstrated its robustness and usefulness for postprocessing the velocity anywhere in the domain. Additionally we have shown that it can be added to an existing drop model [27] to accurately model the movement of drops near corners.
Future work could include extending this model to three dimensions. A boundary integral method for three-dimensional drops in free space has been developed [32,33]. In three dimensions Lipschitz domains admit both corners and edges, so the RCIP algorithm described in Section 4.2 must be further developed. In [15] such an approach is used to compute the capacitance of a cube.