An integral equation method for closely interacting surfactant-covered droplets in wall-confined Stokes flow

A highly accurate method for simulating surfactant-covered droplets in two-dimensional Stokes flow with solid boundaries is presented. The method handles both periodic channel flows of arbitrary shape and stationary solid constrictions. A boundary integral method together with a special quadrature scheme is applied to solve the Stokes equations to high accuracy, also for droplets in close interaction. The problem is considered in a periodic setting and an Ewald decomposition for the Stokeslet and stresslet is derived to make the periodic sums convergent. Computations are sped up using the Spectral Ewald method. The time evolution is handled with a fourth order, adaptive, implicit-explicit time-stepping scheme. The numerical method is tested through several convergence studies and other challenging examples and is shown to handle drops in close proximity both to other drops and solid objects to a high accuracy.


Introduction
The study of deforming droplets on the micro scale is motivated by several applications, one being the design of lab-on-a-chip-devices. In many cases, the study of deformable droplets in a confined flow is especially important. For example, the channel geometry can be used to control the behaviour of the droplets; regarding transport, splitting and fusing of the droplets [7]. It is also of interest for the study of flow through porous media, which is relevant to many industries such as e.g. oil recovery. A review of the physics of the problem is given by Zhang et al. [47].
On this micro scale, the flow can be modelled by the Stokes equations. The surface area to volume ratio is typically very high and interfacial forces are important for the flow dynamics. Surfactants are molecules that alter the surface tension of a drop, which changes the interfacial dynamics and thus the behaviour of the whole system. The inclusion of surfactants is an important tool in drop creation and coalescence prevention [3]. A review of drops and bubbles in shear flow without constrictions is given by Rallison [39]. Regarding the flow through straight capillaries, Olbricht and Kung [32] performed an extensive experimental study of the shape of a drop as a function of several physical parameters, such as Capillary number and viscosity ratio. Furthermore, Shapira and Haber [40] studied how the shape of a droplet between two parallel plates is affected by the ratio of droplet diameter and channel height using small deformation analysis and reflections. The study of physical parameters is clearly of interest, however this paper focuses on the development of a highly accurate numerical method for simulating deforming droplets in wall-bounded Stokes flow with the inclusion of stationary, solid constrictions.
The numerical method described in this paper contains a boundary integral equation method. An overview of some studies of deforming droplets with different kinds of constrictions using boundary integral methods follows: In 3D, Zinchenko and Davis [51] developed a method for simulating one deforming drop squeezing through solid particle constrictions. This method is based on a Hebeker representation for the solid particle contribution and formulates a system of Fredholm integral equations of the second kind. They considered a uniform flow pushing the drop through constrictions and studied the effect of Capillary number and viscosity ratio between drop and bulk on the deformation. To handle interactions between droplets and particles, a special desingularization technique was used. Furthermore, in [11] the authors considered the motion of a single drop between two parallel plane walls, but for this paper the solid walls were handled by modifying the Green's functions. Janssen and Anderson [16] regarded the deformation of drops with unity viscosity ratio between parallel plates, and how the degree of confinement affects the behaviour of the drop. This method is based on a boundary integral formulation in 3D, where the walls are taken into account by using Green's functions associated with the walls. They have since extended their method to handle also non-unity viscosity ratios [17] and also regarded unity viscosity ratios of drops together with insoluble surfactants in [18]. This method is however restricted to flat parallel walls. Tsai and Miksis [44] studied the dynamics of a 3D axisymmetric drop in straight and capillary tubes as a function of viscosity ratio and Capillary number.
In 2D, Zhou and Pozrikidis [48] considered suspensions of drops in channels, using a periodic suspension of viscous drops. Here the drops were ordered in a single file, but studies were also made for random suspensions [49]. A Fredholm integral equation of the second kind was obtained using periodic Green's functions that represented the flow due to a periodic array of 2D point forces in a channel. In the above references, the flow was driven by the relative translation of the two walls. In [50], the flow was driven by a constant pressure drop. Li and Pozrikidis [24] considered wall-bounded channel flow of a suspension of many droplets, for different Capillary numbers and viscosity ratios. They used a boundary integral equation with modified periodic Green's functions to take the walls into account. DeBisschop et al. [8] considered the motion of a two dimensional bubble rising in an inclined channel in Stokes flow. They considered both clean and surfactant-covered bubbles. The fluid velocity was computed using a periodic Green's function (in the x-direction). The authors compared their results with that of experiments regarding inclined walls, showing good agreement.
Other numerical methods to simulate deformable droplets in wall-bounded flow include Lee and Pozrikidis [23], who considered the effect of surfactants on the deformation of drops and bubbles in flow with non-zero Reynolds number. They used finite-differences for the Navier-Stokes equations, finite volume for the insoluble surfactants and Peskin's immersed interface method for the interface tracking. In 3D, Wang and Dimitrakopoulos [46] studied deformable drops in a square channel using a spectral boundary element method. Mortazavi and Tryggvason [29] studied three dimensional deforming drops in a tube using a finite difference/front-tracking scheme. In both 2D and 3D, Claus and Kerfriden [6] used a cut finite element method to simulate deforming bubbles in Navier-Stokes flow. Their results included those of a bubble squeezing through a 5 : 1 : 5 contraction/expansion micro-channel. Chung et al. [5] investigated the effect of viscosity ratio and Capillary number on a similar construction, using the finite element front-tracking method.
Surfactant-covered droplets in two-dimensional Stokes flow without the presence of walls and solid constrictions have been previously simulated with boundary integral equation methods, e.g. by Kropinski and Lushi [22] and Pålsson et al. [34]. The results in [34] were thoroughly validated using exact and semi-analytical solutions obtained by conformal mapping theory. The validation tests showed the ability of the method to obtain highly accurate (e.g. 8 correct digits) in solutions also after a long time with significant droplet deformation and close interactions between droplets.
Resolving the interactions of droplets in close proximity is a challenge for all numerical methods; grid based methods face the need for fine meshes and remeshing, whilst boundary integral equation methods necessitates the handling of nearly-singular integrals. The work in [34] utilised a special quadrature scheme [31] in order to resolve these interactions and achieve very accurate solutions as discussed above. The interactions of other objects in flow, such as vesicles, yield similar challenges. Rahimian et al. [38] used a boundary integral method in 2D extending that by Kropinski [21] to study how vesicles deform over time in confined flows, for example when squeezing through constrictions. Quaife and Biros [36] studied vesicles suspended in a viscous Stokesian fluid, including channel constrictions and other solid geometries. This method was revisited in [37] with the inclusion of an adaptive time-stepping scheme. To handle closely interacting vesicles, interpolation was used. Marple et al. [28] simulated vesicles in periodic channel flows of arbitrary shape, using fast direct solvers. Here, a globally compensated trapezoidal rule was used to handle vesicles in close proximity.
In this paper, the fluid flow problem is considered in a periodic domain. This, however, implies that all periodic images need to be considered when evaluating the fluid velocity. To make this computationally viable, a fast method is needed. For periodic systems, an Ewald summation method [9] is especially suitable. In 3D, Ewald decomposition using different "screening functions" has been performed, see e.g. [4,12,25]. In 2D, Van De Vorst [45] derived the formulation to split the Stokeslet and the stresslet to compute the flow field of a fluid in a domain with pores in the Stokes regime. However, this specific derivation yields a non-symmetric expression for the stresslet decomposition, similar to that obtained by [27] in an alternative derivation for the three dimensional case. The computations of the decomposed expressions can be sped up using the Spectral Ewald method [1,26]. When considering instead a problem without periodicity, either the Spectral Ewald method for freespace [2] or a Fast Multipole Method (FMM) can be utilised. The method in this paper can easily be modified to the free-space case, similar to that in [34].
This paper presents a highly accurate numerical method for the simulation of deforming droplets in a periodic two-dimensional Stokes flow, in the presence of constrictions and channels. The method handles close interactions between drops as well as drops and solid objects without an increase in error. The drops may be clean or covered by insoluble surfactants. The method is an extension of that in [34]. In this paper, the tools described in that paper are extended upon to include also solid walls and stationary objects. The method in this paper is general, i.e. considers both drops and solids in any configuration, and makes no distinction between channel walls and other solid objects. Moreover, this paper also includes a new derivation of the split of both the Stokeslet and the stresslet, using the Hasimoto screening function. This gives a decomposition of the Stokeslet equal to that of Van De Vorst, but a symmetric expression for the split of the stresslet.
The problem setting in this paper is limited to two dimensions. Certain physical effects will be lost through this simplification, however it has been noted that a substantial degree of physical relevance remains [21,48]. Furthermore, the reduction in dimension allows for larger simulations with an increased number of close interactions, due to the substantially reduced computational cost. In terms of numerical methods, this work shows the advantage of boundary integral formulations for highly accurate treatment of interface dynamics and close interactions, and there is ongoing work to develop the same abilities in three dimensions, see e.g. [41] and the references therein.
The paper is organised as follows: in §2 the governing equations, nondimensionalisation and boundary integral formulation are introduced. The complete numerical method is described in §3, and in §4 the Spectral Ewald method in two dimensions to handle periodicity is described, together with the decomposition of the Stokeslet and the stresslet and truncation error estimates to facilitate parameter selection. In §5 the capabilities of the method are demonstrated through numerical tests.

Problem formulation
The equations and mathematical tools needed to simulate surfactant-covered drops in a wallbounded flow with solid constrictions are described in this section. First, the Stokes equations which govern the flow of the problem and the convection-diffusion equation for the surfactant concentration are stated both in dimensional and nondimensional form. Then follows a description of how to reformulate Stokes equations for deformable drops in the presence of stationary solid objects and walls into an integral equation. Finally, the periodic extension of the problem is described.

Governing equations
The equations governing the physical problem are the incompressible Stokes equations, which in their dimensional form are Here, Ω 0 is the bulk fluid surrounding the drops and solid objects and Ω k is the interior of the drop k. There are in total N Γ drops. Furthermore, p k is the pressure and u k the velocity, for k = 0, 1, . . . , N Γ . The drops and the bulk fluid are separated by the interfaces Γ k , on which the normal stress balance holds, wheren k is the outward facing normal, e k the strain tensor for the bulk and the interior of the drops, κ k = ∇ s ·n k the curvature and ∇ s the surface gradient for s traversing the drop k in an anti-clockwise direction. Moreover, σ k := σ k (s, t) is the surface tension coefficient of drop k at time t. The fluid velocity is continuous on the drop boundaries, i.e. u k = u 0 on Γ k . The interfaces are discretised anti-clockwise with a parameter s ∈ [0, L k (t)] where L k (t) is the length of the interface of drop k at time t. An example of a domain configuration can be found in Figure 1. The drops translate and deform according to the ODE for all points x k ∈ Γ k . The boundaries of the solids are denoted by γ k , for all solids k = 1, . . . , N γ . All solid boundaries have a no-slip boundary condition, i.e. the fluid velocity relative to the solid boundaries is always zero. The flow problem to be considered is the case where deformable drops are moving in channels or close to solid objects in an added flow field u ∞ .
In addition, insoluble surfactants are considered. Their concentration is described by ρ k (s, t) and governed by a convection-diffusion equation on each interface; for x k (s, t) ∈ Γ k , where D Dt is the material derivative and D Γ is the diffusion coefficient along the interface [42]. As the surfactants are insoluble, the mass of surfactants is conserved along each interface, for each drop k = 1, . . . , N Γ . The surfactant concentration and the surface tension coefficient are coupled through an equation of state. Here a linear equation of state is considered, for each drop k, where σ 0 is the surface tension coefficient of a clean drop, R is the universal gas constant and T the temperature. This equation of state can be trivially exchanged to others.

Nondimensionalisation
All lengths are nondimensionalised using a characteristic length r 0 , which unless otherwise stated is defined as the radius of the largest drop. The velocity is nondimensionalised by a characteristic velocity U . This velocity is chosen from the imposed far-field flow as U = max(|u ∞ (x)|). Furthermore, the surface tension coefficient is nondimensionalised by the surface tension coefficient of a clean drop, σ 0 . This leads to a characteristic pressure µU r0 and a characteristic time r0 U . Also, the surfactant concentration is nondimensionalised by the initial surfactant concentration on the largest drop of the problem.
For the rest of this paper, all quantities are considered in their nondimensional form. The Stokes equations (1) then read where λ k := µ k µ0 is the viscosity ratio between the fluid of drop k and the bulk. An inviscid bubble corresponds to the limit where λ k = 0. The no-slip condition on the solid boundaries is u k = 0 for all x ∈ γ k , k = 1, . . . , N γ . The stress-balance over each interface Γ k (2) is rewritten as where the Capillary number Ca is defined as Ca = U µ0 σ0 . The convection-diffusion equation governing the surfactant concentration (4) becomes where Pe Γ = r0µ0 σ0DΓ is the Peclet number. Furthermore, (6) becomes where E = RT ρ0 σ0 is the so-called elasticity number.

Boundary integral formulation
A thorough derivation of the formulation for drops and solid particles in 3D can be found in [51]. Here, the same approach is followed but the formulation is rewritten for the 2D case.
For any point x in the bulk fluid Ω 0 , the velocity can be written as where f k := 1 Ca (σ k κ knk − ∇ s σ k ) from (2), S and D stand for the single-layer and double-layer contributions respectively and β(x) stands for the solid-particle contribution as discussed below. The single and double layer potentials are defined as for a drop interface or solid boundary, where Λ = Γ k for a drop k or Λ = γ k for a solid k. Here, G and T are the Stokeslet and stresslet respectively, which in 2D are defined as where r = x−y. The solid-particle contribution can be defined as a single-layer potential over the solid boundaries, however, this generates an ill-conditioned system. In order to obtain a well-conditioned system the same approach as by Zinchenko and Davis [51] is taken, where a Hebeker-representation [13] is used to represent the solid-particle contribution. In this representation, the flow exterior to the solid particles is represented as a combination of single and double layer potentials, where q is the so-called Hebeker density and η is a proportionality factor which is set to η = 1. Another option to achieve a well-conditioned system would be to use a completion flow as in [35]. Finally, taking the limit of (11) as the point x → Γ and x → γ , gives a system of Fredholm integral equations of the second-kind [51], where x ∈ Γ for all fluid interfaces Γ ∈ Γ := NΓ k=1 Γ k and for x ∈ γ := Nγ k=1 γ k (all solid boundaries). This is a system of Fredholm integral equations of the second kind which needs to be solved for each time step to obtain the velocity u with which the drops are moving.

Periodicity
In this paper, the flow problem is considered with periodic boundary conditions in both the x-and y-direction. When computing the flow u at any point x ∈ Ω 0 ∪ γ ∪ Γ through (17) and (18), this means that D Λ [u](x) and S Λ [u](x) contain the integrals over all periodic images over surfaces Λ. With their periodic replicas, they become where τ (p) = (p 1 L 1 , p 2 L 2 ) T for a periodic box of size L 1 × L 2 and p = (p 1 , p 2 ) T , p 1 , p 2 ∈ Z.
In order for the evaluation of these integrals to be computationally possible, a fast method is needed. In this paper the Spectral Ewald method [1,2,26] will be applied, which makes computing u at N target points with O(N log(N )) cost possible, where N is the total number of discretisation points in the system. Details of this method will be covered in §4.

Numerical method
To compute the evolution of deforming surfactant-covered drops, (3) and (9) need to be solved, generating the following system for all x ∈ Γ k , for k = 1, . . . , N Γ . The velocity u k is determined by solving the boundary integral formulation described in §2. Several components are needed to obtain an accurate solution to this system, most of them are described in detail in [34]. In this section, an overview of the method will be given. Moving the drop interfaces by (21) can result in a clustering of discretisation points on the interfaces. This is not ideal, as it necessitates remeshing. One can instead modify the velocity, as was done by Hou et al. [15] for elasticity problems and Kropinski [21] for drops and bubbles. This approach will instead move the drops with velocityũ k , i.e.
where the normal component ofũ k , u n , is the same as for the fluid velocity u k and the tangential velocity is modified as described in §3.4 such thatũ t = u t . Inserting this new velocity into (21) and (22) and expanding the material derivative gives the new system is the length of drop interface k at time t. A hybrid method for discretising the equations in (21) and (22) is used. The hybrid method consists of two discretisations: "Grid 1" is a panel-based composite 16-point Gauss-Legendre discretisation used on both drops and solids, and "Grid 2" is a uniform discretisation in arc-length, used only for the drops. Grid 1 will be used to determine the velocity from the boundary integral formulation in (17), u k (x, σ, t), which is used to compute the normal component u n used to move the drops. Grid 2 will be used for determining the appropriate tangential velocity, for solving the surfactant equation (24) and for updating both ρ and x in time. To go from the equidistant discretisation to the Gauss-Legendre one nuFFT Polynomial interp. a non-uniform FFT is used, see Greengard and Lee [10]. To go the opposite way, a 16-point polynomial interpolation on each panel is used. A schematic is shown in Figure 3.
The problem setting is as follows; there are N Γ drops and N γ solids. A solid k is discretised by M γ k Gauss-Legendre points (Grid 1), giving a total of discretisation points on the solids as This makes the total number of discretisation points of the system M = M γ + M Γ , and the number of unknowns 2M .

Complex variable notation
When regarding the problem of deforming drops in 2D, it is beneficial to consider the formulation in complex variable notation where a point x corresponds to z = x + iy. Considering z, τ ∈ C, the complex counterparts of S Λ [g](x) and D Λ [g](x) are denoted S Λ [g](z) and D Λ [g](z), for Λ = Γ (layer potential over drop interfaces) and Λ = γ (layer potential over solid boundaries). Furthermore, in complex notation and For the remainder of this section, the problem will be treated in this complex setting.

Computing the velocity u(z, σ, t)
To compute the velocity the composite Gauss-Legendre discretisations ("Grid 1"), of drop and solid interfaces are used. The system to discretise is given by (17) together with (18), using the periodic expressions for S and D in (19) and (20). This system of Fredholm integral equations of the second kind will be solved by a Nyström method. In its discretised version this becomes, for u ≈ u(z ) and q ≈ q(z ), where z are the Gauss-Legendre discretisation points with associated weights w , on either a drop interface or solid boundary, Here, f : (25) respectively, including the periodic extension as mentioned in §2.4. In §4 the computation of these periodic sums is described. Define τ (p) = p 1 L 1 + ip 2 L 2 for p = (p 1 , p 2 ) with p 1 , p 2 ∈ Z and L 1 , L 2 the sides of the periodic box. The periodic double and single layer potentials are then defined as for Λ = Γ, γ and where (for z = z m ) Observe that the limits when z m = z are finite for M (1) , M (2) and M (3) and given by The second integral of (25) corresponding to the integral of M m is more complicated as the kernel is non-smooth. Also, note that special quadrature is needed for the nearly-singular case when drops get close to each other or solid objects/walls. How to deal with both these issues is described in §3.3. In this subsection, it is assumed these problems can be dealt with efficiently and a highly accurate solution obtained for all points on all drops.
The discretised system in (27) is then solved using gmres. The periodic sums are all computed with the Spectral Ewald method, as described in §4. The discretised system in (27) has a unique solution by the Fredholm Alternative, and as it is a Fredholm integral equation of the second kind has spectral properties which enable gmres to converge in few iterations. The authors observe that the number of iterations of gmres vary with viscosity ratio, but have in none of the test simulations felt the need for a preconditioner. 1,4], in §3.2 become near-singular, i.e. when evaluating them at a point z such that z m − z 1, for some m, yielding large numerical errors when the regular Gauss-Legendre quadrature rule is applied. This is the case for example when drops get close to each other or solids. How these errors behave and can be estimated was studied in [34]. To obtain accurate approximations of the integrals at any distance from the interfaces, a special quadrature scheme will be employed. The main idea of the specialised quadrature has been described in [30,31,34]. This quadrature scheme will be employed for near-interactions in the case of integrals of M (k) m for k ∈ [1,4]. In the case of M (4) m special treatment is needed also for the on-surface evaluations. For all kernels in (30), the main idea of the special quadrature is similar. In short, for an integral of the form Γ f (τ )K(τ, z)dτ , the idea is to express the function f (τ ) as a polynomial in τ , where the coefficients can be computed using a Vandermonde system. One can then use recursive formulas to compute the integrals as needed analytically. Following the notation in [30], all the integrals of M (k) m , k ∈ [1,4] can be written on one of the following forms

Special quadrature
for any smooth boundary Λ. How to handle I 1 and I 2 is described in [31,34]. The special quadrature to deal with I 3 can be found in [14,30]. A brief overview of all three cases is given here. Note that the special quadrature treatment is strictly short-ranged, so there is no need to involve any periodicity in the calculations, except when considering drops and solids close to the edge of the periodic box. Furthermore, note that the third integral I 3 can be rewritten in the following way wheren τ is the normal of Λ at point τ . For all three integrals, the approach is the same. Consider Λ as a panel on either Γ or γ, rotated and scaled such that its endpoints are at −1 and 1 in the complex plane. The evaluation point z is rotated and scaled along with Λ. Now, expanding h(τ ) as a monomial, coefficients c can be computed such that where n is the number of Gauss-Legendre points on each panel, here set to n = 16. Inserting the interpolation into the integrals I 1 , I 2 and I 3 the following expressions are obtained where The analytical integrals p , q and r can be computed through recursion formulas, given in Appendix C. The polynomial coefficients d in (32) are given by a polynomial expansion of the function f (τ ) = h(τ )/n τ , similar to that of c .
Regarding the on-surface evaluation of I 3 (z), the same approach is used. It can be seen as computing the quadrature weights for the log-kernel with a particular target point z. Using this approach, the quadrature weights can be precomputed and saved for all quadrature points on a panel in a n × n matrix. This can be extended to a rectangular matrix to include special quadrature treatment also for target points on the neighbouring panels.

Modifying the tangential velocity
As previously stated, the tangential velocity needs to be modified in order to avoid clustering of discretisation points on the interfaces. Denote the velocity computed in §3.2 by u(z) and the modified velocity byũ(z). Both velocities are considered on the uniform grid ("Grid 2"). The normal components of the velocities are the same, i.e.
The tangential component is modified according to (see [31]) Note that the modified velocityũ(s) = [u n (s) + iũ t (s)]n(s) still fulfills the kinematic condition.

Solving the surfactant equations
The equation for insoluble surfactants is solved as described in [34]. The equation to compute the surfactant concentration is given by (24). This equation can be solved using a pseudo-spectral method which generates a system of ODEs to solve, one for each Fourier coefficient of ρ k (α, t); ρ k j (t): for all drops k. Here, f I corresponds to the part of (24) that needs to be treated implicitly due to stiffness, and f E corresponds to everything else which is handled explicitly, i.e.
Also, f k E j and f k I j correspond to the jth Fourier coefficient of f k E and f k I respectively. Since this is a pseudo-spectral method, the surfactant concentration is computed on the uniform grid ("Grid 2").

Spatial adaptivity
As the drop deforms over time the uniform discretisation will keep its points equidistant with distance ∆s due to the modified tangential velocity in (34). As the interfaces can stretch and contract over time, it is beneficial to allow for spatial adaptivity keeping ∆s constant. Using FFTs, this operation is trivial on the uniform grid.

Time-stepping scheme
The coupled system to time-step is where f k E j should be treated explicitly and f k I j implicitly. In [34], a second order time-stepping scheme was utilised. This time-stepper was chosen after a comparison of several schemes in [33]. It has for this paper been updated to a fourth order time-stepping scheme, which gives a considerable gain in computational cost. The scheme needs to handle adaptivity in time for both surfactant concentration ρ and position z, and utilise the same stages for both equations. Therefore, a fourth order adaptive scheme by Kennedy and Carpenter [19] is used, which uses the "ARK4(3)6L[2]SA-ERK" for the explicit parts together with the diagonally implicit "ARK4(3)6L[2]SA-ESDIRK" for the implicit part. These are additive Runge-Kutta methods where adaptivity is acquired by comparing to a lower order scheme. For a Butcher tableau of the scheme, the reader is referred to Appendix C in [19]. The time-step is modified using where s f = 0.8 is a safety factor, tol the given time-stepping tolerance, is machine epsilon and r = max(r z , r ρ ), where r z is the measured error in z and r ρ the measured error in ρ.

Summary of the numerical method
Above, each step of the numerical method is described. Here follows an overview of how the different parts are put together.
Initially, all drop boundaries are discretised with a discretisation that is uniform in arc-length ("Grid 2"). The solid boundaries are discretised with a composite 16-point Gauss-Legendre scheme ("Grid 1"). Surfactant concentration is initialised on the uniform grid of the drop interfaces. Timestepping is performed as described in §3.7. For every stage fromt tot + c dt in the time-stepping scheme, the following steps are taken: 1. The uniform drop discretisation is transformed to the panel-based G-L quadrature through a nuFFT, i.e. "Grid 2" → "Grid 1". 2. The velocity u for timet is computed by solving the integral equation on both the drop interfaces and the solid boundaries, see §3.2. To compute the velocities for the periodic problem Spectral Ewald as described in §4 is utilised. Special quadrature as in §3.3 is used to obtain high accuracy for all discretisation points. 3. Once the velocity u for timet is obtained, this is interpolated back to the uniform grid, "Grid 1" → "Grid 2", for the drop discretisation points. 4. The velocity is then modified toũ using (34). Additionally, f E (t) as in §3.5 is computed. 5. The new position and surfactant concentration at timet + c dt is computed using the method in §3.7. 6. The surface tension coefficient at timet + c dt is computed through (10). 7. If an interface length has changed sufficiently, the number of discretisation points is modified through FFTs to keep ∆s constant, see §3.6. Also, a Krasny filter of level 10 −12 is applied to both position and surfactant concentration.

Periodicity and Spectral Ewald
As mentioned in §2.4 the flow problem is considered in a periodic setting in both x-and y-direction. The integrals to compute are (19) and (20) for the Stokeslet and stresslet respectively. They are discretised using the Gauss-Legendre quadrature described in §3, with quadrature nodes and weights x n , w n , n = 1, . . . , M Λ , for M Λ the total number of discretisation points on a boundary, Λ. Thus, the approximations of the integrals (19) and (20) are As previously in §2.4, τ (p) = (p 1 L 1 , p 2 L 2 ) T for a periodic box of size L 1 × L 2 and p = (p 1 , p 2 ) T . The Stokeslet G jl and the stresslet T jlm are defined as in (14) and (15) respectively. They both decay very slowly and the sums (37) and (38) are only conditionally convergent. Ewald decomposition [9] is used to remedy this, where each sum is split into two parts: one which contains the singularity and converges rapidly, referred to as the "real space" sum, and one which contains a smooth periodic function, and thus converges quickly in Fourier space; the "k-space" sum. In order to compute the "k-space" sum at

Ewald decomposition
To illustrate the idea of Ewald decomposition, the split into "real space"and "k-space" is first computed for the Green's function for the biharmonic equation. In 2D, this Green's function has the following form and it is the fundamental solution to −∆ 2 B(|x − y|) = δ(|x − y|). The choice of constant α is free, and is here chosen to α = 3 2 preserve the 3D relation between the Green's function and the Stokeslet, i.e. G jl (r) = (∆δ jl − ∇ j ∇ l ) B(|r|) [2].
When considering a sum where the asterisk in the first sum corresponds to the exclusion of the term x − x n − τ (p) = 0, the aim is to find a split into B R (r, ξ) and B F (k, ξ) such that Here, the first sum corresponds to the "real space" sum, and the second to the "k-space" sum. Since the term where x − x n − τ (p) = 0 is excluded from the sum u B , it should also be excluded from the "k-space" sum when evaluating at a target point x = x n , where x n is any of the source points. This can be done by adding the limit to the expression above. .
The screening function should be defined such that B R (|r|, ξ) is short-range and B F (|r|, ξ) is smooth and long-range. The Hasimoto screening function, defined as meets these criteria, where r = |r| and k = |k|. In Fourier space, the convolution B(|r|) * γ(r, ξ) is computed as a multiplication, which gives B F (|k|, ξ) = γ(k, ξ) B(|k|, ξ). For the biharmonic Green's function, B(|k|) = −1/k 4 . How to compute B R (|r|, ξ) by convolution is described in Appendix A.1. The split obtained is the following,

Stokeslet
To find the split of the Stokeslet, first note that the Stokeslet can be expressed as an operator K jl acting on B(|r|) [2], where K jl = ∆δ jl − ∇ j ∇ l . Thus, to compute the "real space" part of the Stokeslet, G R jl , this operator is applied to B R which gives wherer j = r j /r. Similarly, the Fourier space part can be expressed as G F jl (k, ξ) = K jl B F (|k|, ξ), where K jl denotes the pre-factor that is produced when K jl is applied to e −ik·r , i.e. K jl = −δ jl k 2 + k j k l . This gives wherek j = k j /k and k = |k| for k = (k 1 , k 2 ). Thus, u G j (x) in (37) can be written as In the case of the target point x = x n for any source point x n , the self-contribution that arises from the Fourier part needs to be removed. Thus, needs to be added to the expression, where γ is the Euler-Mascheroni constant.

Stresslet
The Ewald decomposition of the stresslet is computed in a similar manner to the Stokeslet, by applying an operator K to the decomposition of the biharmonic Green's function. The operator relating the stresslet and B(r) is [2], Thus, the "real space" part of the stresslet is given by By applying K jlm = −i (δ jl k m + δ jm k l + δ lm k j ) k 2 − 2k j k l k m to B F (|k|, ξ), the "k-space" part of the stresslet is computed as The complete expression thus reads For the stresslet, there is no "self-interaction term" as is the case for the Stokeslet, as lim |r|→0 T R jlm (r, ξ) − T jlm (r) = 0. The term T F,(0) is chosen to guarantee zero-mean flow through the primary periodic cell, and also ensures that the stresslet identity is met [1]. This corresponds to setting T F,(0) jlm (y) = 4πδ lm y j .
As a side note regarding the stresslet, it can also be computed using derivatives of the Laplace Green's function, L(r), and the Stokeslet, i.e.
Computing T R (r, ξ) and T F (k, ξ) through this relation using the Hasimoto screening of L(|r|) into L R (|r|, ξ) and L F (|k|, ξ) and G R jl (r, ξ), G F jl (k, ξ) from above, generates the same results as those in (45), (46). The split of L is derived in Appendix A.2. Moreover, if instead using the so-called Ewald screening function to split L, while keeping the Hasimoto screening function for the Stokeslet, G, this generates the same expression as that found by Van De Vorst [45]. This expression, however, is not symmetric, and will not be used in this work.

Truncation errors
The sums in (44) and (48) converge fast, but cannot be computed numerically without truncation. To decide where to make the truncations, the "real space" sum and the "k-space" sum need to be regarded separately and their truncation errors estimated. The error estimates and their derivations are described in detail in Appendix B and inspired by the error estimates in [2]. Here, the case of a square box, i.e. L 1 = L 2 = L is considered for simplicity.

Truncation errors for the "real space" sum
The two "real space" sums for the Stokeslet and stresslet respectively are defined as For these sums, only target points within a cut-off radius defined as r c will be considered. Comparing u G,R j and u T,R j with the truncated sumsũ G,R j andũ T,R j , the RMS of the truncation error is given by for the Stokeslet and stresslet respectively. The truncation errors are estimated as for the Stokeslet, and for the stresslet, where The derivations of these estimates can be found in Appendix B.1. The truncation errors and estimates for both Stokeslet and stresslet are shown in Figure 3. Comparing the "real space" estimates to those empirically obtained for 3D in [1,25], the estimates show the same asymptotic behaviour, but vary in the constant factor in front and powers of (ξr c ) needed.   (50) and (51). The system is Ns = 10 3 randomly distributed point sources within a square of size L = 2π, and Nt = 10 2 randomly distributed target points in the same square.

Truncation errors for the "k-space" sum
For the "k-space", the sums to compute are defined as and they are truncated in "k-space"such that k = (k 1 , k 2 ) for k 1 , k 2 ∈ [−k ∞ , k ∞ ]. The RMS of the truncation error is estimated as for the Stokeslet and for the stresslet, with Q G and Q T as defined previously. These estimates are derived in Appendix B.2. The truncation errors together with the estimates are shown in Figure 4 for both the Stokeslet and the stresslet. Additionally, comparing the "k-space" estimates to those for 3D in [1,25], the estimates again show the same asymptotic behaviour, but the expressions in front of e −k 2 ∞ /4ξ 2 differ. The "k-space" truncation estimates are not as precise as their "real space" counterparts, but always overestimate the errors.

The Spectral Ewald method
Using Ewald decompositions such as those in (44) and (48), the sums to compute are now rapidly converging. For a system with N discretisation points, computing the sums directly results in an O(N 2 ) complexity. To speed this up, the Spectral Ewald method [1,2,25] is used, which makes the computations of the k-space sums O(N log(N )) in cost. The method is thoroughly described in the previous references for the three dimensional case, and generalises to 2D easily. In short, the key approach is to evaluate the Fourier space sums on a grid in the "k-space" sum which enables the use of FFTs of the size M 2 to speed up the computations. To spread source and target points to the grid truncated Gaussians with P 2 points support are used, and the shape parameter of the Gaussians is determined to minimise the approximation error for this given P . The real space sums can be computed in O(N ) time, by constricting the evaluation only to points in a near-neighbour list of each point x t , defined as Under the assumption of a constant number of those near neighbours, to create such a list is also O(N ). The parameter ξ from the screening function, decide how much work is put into the "real space" sum and "k-space" sum respectively.
There are several parameters to set in the method. To keep the number of nearest neighbours constant in the real-space sum as the system is scaled up, the cut-off radius r c is set first. Using the estimates in (50) (Stokeslet) and (51) (stresslet), for a given tolerance tol e the splitting parameter ξ is computed. From ξ the corresponding k ∞ = M 2 is computed from the estimates in (52) (Stokeslet) and (53) (stresslet). Moreover, P is set to P = 24. This keeps the approximation errors of the method close to round off and is in 2D not very costly.

Numerical results
The numerical method concerning the drop deformation has previously been thoroughly validated using conformal mapping techniques in [34]. In this section the extended numerical method described in §3 including solid objects and walls is thoroughly tested, through convergence studies and difficult test cases. Each case is described in detail below.

Convergence study -drop squeezing through a constriction
Firstly, a convergence study of drops of different viscosity ratios squeezing through a solid constriction is performed. The domain at time t = 0 consists of two stationary solid discs of radius r = 1 that are placed with centre points in (±1.375, 0), i.e. their minimum distance is = 0.75r. Also, a drop with radius r and viscosity ratio λ is placed with its centre point in h = (0, 2.1). The drop is pushed down through the constriction by a imposed flow, u ∞ = (0, −1). The Capillary number is set to Ca = 1. The viscosity ratios investigated here are λ = 0.5 and λ = 2. The domain set-up is shown in Figure 5.  Comparing the solutions. Denoting the reference solution by z ref , and a coarser solution using N discretisation points for τ , the aim is to compute the difference between z ref and τ . Previously (see e.g. [34]), the coarser solution has been upsampled to the same size as z ref using FFTs. For this convergence study, however, this will not be an optimal approach as will be explained in §5.1.1. Here, instead a different approach based on a normal projection onto the reference solution will be considered. In essence, the closest distance between a point in the coarse discretisation and the reference solution needs to be measured. A potential tangential shift of the point is irrelevant. A schematic of this procedure can be found in Figure 6. Considering the N ref discretisation points of z ref , they are represented on the equidistant grid and can be seen as discrete points of a periodic function z(α), where α ∈ [0, 2π]. It is therefore possible to obtain their Fourier coefficients through an FFT. Once these coefficients have been obtained, a normal projection of a coarse discretisation point τ k onto the reference interface can be found through a minimisation procedure. This procedure can be formulated as finding theα such that τ k − z(α) ∞ is minimised. The difference between the reference solution and the coarse discretisation point τ k is defined as this distance. To compute the difference between a coarse solution and the reference solution, this procedure is repeated for all discretisation points of the coarse solution.
z ref τ α z(α) τ k Figure 6: Schematic of how the difference between reference solution and coarse solution is computed. Normal projection z(α) of τ k from coarse discretisation onto the reference interface is marked in red.

Clean drops
In the case of no surfactants, the evolution of the domains and the error as a function of 1/∆s can be seen in Figure 7 for λ = 0.5 and in Figure 8 for λ = 2. The simulations were run to two timestepping tolerances tol 1 = 10 −6 (marked with red dashes) and tol 2 = 10 −8 (marked with red dots). Throughout the simulations ∆s is approximately constant in time, due to the spatial adaptivity. Each time t ∈ [1, 2, 3, 4, 5] is represented by a black line, with diamond markers for tol 1 and square markers for tol 2 , showing the error as a function of 1/∆s at that particular time t. A comparison between the two viscosity ratios shows that for the same non-dimensional instance in time, the lower viscosity drop has deformed more, which is to be expected. For both cases the error decreases with an increase in 1/∆s until the time-stepping error dominates. The error is roughly the same for all times t. For both λ = 0.5 and λ = 2, the set tolerance is reached at an approximate ∆s ≈ 0.026 for tol 1 = 10 −6 and ∆s ≈ 0.020 for tol 2 = 10 −8 .  : Clean drop with viscosity ratio λ = 0.5 squeezing through constriction. Left: Grey represents stationary solid discs, blue represents drop at times t = 0, 1, 2, 3, 4, 5. Right: relative error measured in max-norm as a function of 1/∆s for two time-stepping tolerances: tol 1 is marked with a red, dashed line, and corresponding errors black lines with diamonds (♦), tol 2 is marked with a red, dotted line, and corresponding errors black lines with squares ( ).
The influence of time-stepping tolerance. It is clear from Figure 7 and Figure 8 that both the correct tolerances can be reached. One would expect the less accurate solution to be cheaper to compute, however this is generally not the case. In Table 1 chosen such that both tolerances are reached, and correspond to 1/∆s ≈ 76 and 1/∆s ≈ 204. The cost of a simulation is defined as the number of velocity computations, i.e. the number of integral equation solves, since this is the most expensive part of the algorithm. From Table 1 it is clear that the difference in cost between the two set tolerances is negligible, for both values of ∆s. Moreover, the cost for the larger tolerance tol 1 is even slightly higher than for tol 2 , in the case of the smaller ∆s. This is explained by Figure 9, where the magnitude of each successful time step is shown over time, for both tolerances and ∆s. For larger values of ∆s, i.e. coarser discretisations, the larger tolerance tol 1 allows for larger time steps to be taken. However, as ∆s is decreased, the time steps become of equal size, except for at the very beginning. The reason for this is discussed in the following paragraph, but one can conclude that using a stricter tolerance infers practically no additional cost.  Table 1: Number of failed and successful time steps up to time t = 5 for a clean drop squeezing through a constriction, for λ = 0.5, using the method and time-stepping scheme described in §3.
The fourth order time stepping scheme described in §3 allows the method to take much larger time steps than the previous second order method in [33]. See Table 2 for an overview of the number of time steps taken with tol 1 and 1/∆s ≈ 76. The number of time steps taken increases from 191 with the fourth order method to 3685 with the second order method. This corresponds to approximately six times as many velocity evaluations. However, the required time steps are larger with the higher order method and they can come close to the stability limit. This is noticeable when regarding the equidistant spacing of the discretisation points. When taking very small time steps, such as is the case with the second order method, the points are held equidistant through time. With the larger time steps in the fourth order method, this only holds up to the time-stepping tolerance. This is the reason why computing the errors using FFTs for the coarser solutions and zero-padding is not viable, as it introduces additional errors. How the spectrum looks for the two tolerances is shown in Figure 10. Thus, there is little computational gain when relaxing the time-stepping tolerance.
It is the recommendation of the authors, to in light of this information always run the simulations to the stricter time-stepping tolerance tol 2 . This keeps the discretisation points equidistant with a clean Fourier spectrum (no ringing), and infers practically no additional cost.  Table 1.

∆s
Tolerance #failed dt #successful dt #velocity computations 0.0131 10 −6 1 3685 7370 Table 2: Number of failed and successful time steps up to time t = 5 for a clean drop squeezing through a constriction, for λ = 0.5, using the second order time-stepping scheme described in [33].

Surfactant-covered drops
In this section, the same set-up as above is used with the addition of insoluble surfactants on the drop interface. The non-dimensional initial surfactant concentration is ρ 0 = 1. Furthermore, the elasticity number is set to E = 0.2 and the Péclet number is set to P e = 10. How the drop squeezes through the constriction is shown in Figure 11, together with the evolution of surfactant concentration on the interface. It is clear that the surfactant concentration affect the drop deformation, especially in places of high curvature.
The relative error in max-norm compared to the reference solution is shown in Figure 12. Several things should be noted with these errors. Firstly, the drop and surfactant errors are on different levels. I.e. for a set time-stepping tolerance of 10 −8 , the errors in position will be stable at around 10 −9 whilst the surfactant concentration error level out at approximately 5 · 10 −8 . This could be easily controlled by using different time-stepping tolerances for the two quantities. Secondly, it is also clear that the error is larger for a set ∆s for the times t = 4 and t = 5 than for the earlier times, for both position and concentration. This is due to the increase in curvature of the drop shape, which can be seen in Figure 11 for time t = 4. This is a consequence of the fact that a smaller ∆s is needed to resolve interfaces with high curvature. Practically, this can be handled in a simulation by performing spatial adaptivity not only to keep ∆s constant, but also to decrease it as the curvature increases. This is currently not performed in the simulations.

Multiple drops in a channel
Here, the simulation of multiple drops in a periodic channel is shown. The set-up consists of 15 drops of viscosity ratio λ = 5: two with radius 0.5, six with radius 0.25 and seven with radius 0.15. The walls are parametrised with C ∞ curves and constructed through a superposition of sinus curves. Furthermore, a solid disc of radius 0.5 is placed in the channel. The periodic length is L = 2π. An added Poiseuille flow is driving the movement of the drops. The initial set-up is shown in Figure 13. The minimum distance between the channel walls is 0.45. The evolution of drops from time t = 0 to t = 200 is shown in Figure 14. The drops are initially discretised with 20, 10 and 6 panels for the three different drop sizes respectively, giving ∆s = 0.01. The solids are discretised with a similar ∆s. Through the simulation the Capillary number is set to Ca = 5. During the whole simulation (time t = 0 to t = 200), the minimum distance between two drops is 0.005 and between a drop and a solid 0.008. The time-stepping tolerance is set to 10 −8 and the area error is less than 2.5 · 10 −5 for all times, and can be seen in Figure 15. The increase in area error at time t ≈ 60, is due to the increase in curvature in the yellow drop as seen in Figure 14(c). This higher curvature is due to the large Capillary number chosen for this simulation, and as a consequence more discretisation points to maintain a low error are needed. With twice as many points, i.e. ∆s = 0.005, the area error stays under 10 −7 at all times, see Figure 15. With the addition of surfactants, the surface tension of the drops is lowered and the drops therefore deform more. In Figure 16 the deformation of the surfactant-covered drops can be seen for the case when the simulation in Figure 13 has been modified to include an initial surfactant concentration on all drops ρ 0 = 1, with elasticity number E = 0.5 and Péclet number P e = 1000. A comparison of the deformation for drops with and without surfactants is shown in Figure 17. As can be seen, the addition of surfactants allows the drops to deform more. An example of the surfactant concentration on one drop (the drop marked in yellow in Figure 16) can be seen in Figure 18 (left) for times t = 0, 5, 15, 25, 35. The minimum distance between drops is 0.04 and between drops and solids 0.03. The conservation error of the surfactants and area error of the drops can be seen in Figure 19. It is clear that the surfactants suffer from errors greater than that of the drops position, which was already noted for the previous test case.

Conclusions
An accurate method for simulating droplets together with walls and solid stationary objects in two dimensional Stokes flow has been presented. The method allows for highly accurate solutions due to the boundary integral formulation together with the special quadrature scheme that allows for nearinteraction to be well resolved. The method can handle both channel walls and solid constrictions for flow problems in a two-dimensional periodic setting. To match the high order accuracy in space, a fourth order adaptive time-stepping scheme is utilised.
In order to compute the periodic expressions efficiently, Ewald decompositions for both the Stokeslet and the stresslet have been derived and their truncation errors estimated. The decomposed expressions are then computed efficiently with the Spectral Ewald method. The accuracy of the method has been demonstrated through convergence tests both for clean and surfactant-covered drops. It's stability and robustness have been tested through challenging examples.

Acknowledgements
We are very grateful to Dr. Rikard Ojala for his contributions at the initial stages of this work. This work is supported by the Göran Gustafsson Foundation for Research in Nature and Medicine. A.K.T also gratefully acknowledges the support from the Swedish Research Council, Grant no 2015-04998.

Appendix A. Ewald decompositions
Appendix A. 1

. Decomposition of B(|r|)
In order to split the Green's function for the biharmonic equation, B(|r|) as defined in (39), the following quantities need to be computed: B F (|k|, ξ) and B R (r, ξ). The first correspond to the "kspace" and can be easily computed by where k = (k 1 , k 2 ), k = |k|, γ(k, ξ) is the Fourier transform of the Hasimoto screening function as defined in (40) and B(|k|) = −1/k 4 . The "real space" part, B R , is obtained through for the Hasimoto screening function γ(r, ξ), where r = |r|. Note that B R is expected radial, but no assumption of this is made. Using that it can be written To compute this integral, first switch to polar coordinates (κ, θ), where (k 1 , k 2 ) = κ (cos(θ), sin(θ)) , r = r (cos(ψ), sin(ψ)) .     The integral to compute B R can then be rewritten as which integrated over θ becomes where J 0 (x) is the Bessel function of the first kind of order 0. This integral is difficult to compute, and is therefore differentiated w.r.t. r according to a trick as described in [43]. Note that ∂J0(κr) ∂r = −κJ 1 (κr). Differentiating B R w.r.t. to r then becomes To compute this integral, consider it in two steps: Integrating this w.r.t. r gives which indeed is radial.
Appendix A.2. Decomposition of L(|r|) In 2D, the Laplace Green's function is defined as which is the fundamental solution to −∆L(|r|) = δ(|r|). When considering a sum where the asterisk in the first sum corresponds to the exclusion of the term x − x k − τ (p) = 0, the aim is to find a split into L R (r, ξ) and L F (|k|, ξ) such that where the last term is only included if x = x n for any n ∈ [1, M Λ ]. The Hasimoto split of L is obtained by convolving L with the Hasimoto screening function as defined in (40). This gives that L R (r, ξ) = L(|r|) − L(|r|) * γ(|r|, ξ) L F (r, ξ) = L(|r|) * γ(|r|, ξ).
Using that L(k) = 1/k 2 , it follows that the "k-space" part corresponds to Similarly as for B R in Appendix A.1, L R can be written as The inverse Fourier transform of this is Now, switching to polar coordinates (A.2) and integrating over θ, this reads using the same trick as in (A.3).

Appendix B. Truncation errors
For both the Stokeslet and the stresslet, there will be truncation errors when the infinite sums are truncated for computation. The "real space" sum is only evaluated for point pairs at a distance smaller than some cut-off radius r c , and the "k-space" sum is evaluated for all k 1 , k 2 ∈ [−k ∞ , k ∞ ] for some k ∞ . In general, denoting the truncated computation of one sumũ(x) and the full solution u(x), the RMS truncation error is defined as Appendix B.1. Real space sum The real space sums to compute are for the Stokeslet and stresslet respectively. The truncated sums are denoted byũ G,R j andũ T,R j .
Stokeslet. Starting with the Stokeslet, note that the error of computing such an infinite sum is given by where F L(x) = {(x s , p) : |x − x s − τ (p)| > r c } for a chosen cut-off radius r c . The RMS error is given by |u G,R (x n ) −ũ G,R (x n )| 2 .
With the approximation that E 1 (x) ≈ e −x x for large x, this can be computed and simplified as Again, using a series expansion for E 1 (x) and approximating E 1 (x) ≈ e −x x for large x, δu T,R can be simplified into δu T,R 2 ≈ 2πQ T L 2 ξ 2 r 2 c e −2ξ 2 r 2 c , (B.2) and the error and estimate is shown in Figure 3 (right). Here, for the Stokeslet. All the points x n are contained within a periodic box of size L × L, which means that k ∞ = 2π L M 2 when covering the box with M 2 points in a square grid. The RMS of the truncation error is computed as whereṼ is a circle with radius L/2 containing all source and target points. Corresponding expressions hold for the stresslet.
To estimate the Fourier space truncation error for the Stokeslet, let where the same approximation of the integral as in [20] is used. With G F jl = δ jl −k jkl 1 + k 2 where polar coordinates have been used in the last step. Integrating with respect to θ gives that e jl (r) ≈ π L 2 κ>k∞ 1 + κ 2 4ξ 2 e −κ 2 /4ξ 2 κ 2 J 0 (κr)κdκ, for r = x 2 + y 2 . Using J 0 (x) ≈ √ 2 √ πx for large x, the integral above is approximated as where Γ(ν, x) is the incomplete Gamma function, and Γ( This gives that δu G,F 2 ≈ 4Q G π L 5 k ∞ e −2k∞/4ξ 2 . (B.4) The truncation error and estimate can be seen in Figure 4 (left).
Stresslet. The approach to derive a truncation error estimate for the Fourier space sum of the stresslet is similar to that of the Stokeslet above. First, let u T,F (x) −ũ T,F (x) j = e jlm f lnm where e jlm (r) = 1 L 2 k, |k|>k∞ T F (k, ξ)e ik·r ≈ 1 L 2 |k|>k∞ T F (k, ξ)e ik·r dk.
UsingT F jlm = 1 + k 2 4ξ 2 δ jlkm + δ jmkl + δ lmkj − 2k jklkm e −k 2 /4ξ 2 k , and similar steps as above, with 1 8 2 j,l,m=1 e jlm can be approximated as with the use of E −1/4 (x) ≈ e −x /x for large x. The RMS error then becomes In this last step the same simplifications as for the Stokeslet have been applied, as well as erfc(x) ≈ e −x 2 √ πx for large x. The truncation error and estimate can be seen in Figure 4 (right).

Appendix C. Special quadrature
Here the recursion formulas for computing p , q and r in (33)  Note that if z is within the contour created by the transformed panel Λ and the real axis from −1 to 1, then a residue of 2πi must be added or subtracted from p 0 depending on if z has a positive or negative imaginary part respectively. Similarly, to compute q l , the recursion is q = zq −1 + p , = 1, . . . , n − 1.
Finally, the recursion for r is given by for γ = (τ (2) − τ (1) )/2 where τ (1) and τ (2) are the endpoints of the panel Λ before it was scaled and rotated. This is needed as the log-kernel is not scale and rotation invariant.