Conclusion

γt

Method C, flat triangles Method S, subdivision surface Pozrikidis 2001 [46], curved triangles Le 2009 [50], curved triangles Le 2010 [205], subdivision surface Huang 2012 [160], subdivision surface Zhu 2015 [51], spherical harmonics

(a)G= 0.05,ˆκB= 0.0375.

0 0.05 0.1 0.15 0.2 0.25

0 2 4 6 8 10

G= 0.05 G= 0.2

D

˙ γt

Method A Tsubota [187], model KN Method C Tsubota [187], models H & J

(b)ˆκB= 2/15forG= 0.05andG= 0.2.

Fig. 4.5: Taylor deformation parameter computed via BIM in comparison to data from the literature. The result by our code is represented by Method A and C (the latter matches with B, D and E, compare figure4.4). Reprinted from publication [Pub1] with permission from Elsevier.

4.3.3.3 Conclusion

In the end, the most practical and flexible algorithms for flat triangles appear to be Methods E and C. The first exhibits smaller errors while the latter was found to be more stable in dynamic simulations (which is the reason why we use it in our research projects in chapter6). Compared to Method C, Method B gives slight worse results. Method D appears to be the most robust algorithm (i.e. it is only little affected by the homogeneity of the mesh) and has the best convergence properties, but unfortunately it reaches error levels comparable to the other methods only at very high and thus unpractical resolutions. Method A should not be used to represent the Helfrich model. Still, it might be sensible for stabilization purposes [207]. If the rest of the code (such as the computation of the boundary integrals) employ subdivision surface techniques, Method S can be used without much additional overhead. Hence, different algorithms work better in some contexts while others are preferable in other situations.

5 The boundary integral method:

Algorithmic and numerical treatment of the boundary integral equation

5.1 Putting the boundary integral method into perspective

Other methods Boundary integral methods (BIM) are one of the oldest and most popular methods to study the dynamics of suspended objects in viscous flows [28,42, 131,134], starting with an article by Youngren and Acrivos [135] in 1975. But of course, it is not the only one.

In recent years, the Lattice-Boltzmann [208–210] in conjunction with the immersed boundary method [26,34,77,78,167,211–213] has also become an often used method owing to the comparably simple implementation, incorporation of inertia and the capability of straightforward large-scale parallelization. Alternative mesoscale approaches such as multiparticle collision dynamics [214–219] (also called stochastic rotation dynamics) and (smoothed [11,61,220– 222]) dissipative particle dynamics [33,94,182,185,223–226] are also often and actively used in the biofluidic community. Other methods include the finite volume method [61,99,227–229], the immersed finite element method [230, 231], the moving particle semi-implicit method [232]

and volume of fluid [233]. Also see recent reviews [28,65,234,235].

Advantages of BIM Compared to other algorithms, the boundary integral method has various advantages and disadvantages. One of the most important general advantages is the fact that only the surfaces of the objects occur in the method. This means that the resolution of the surfaces is not tied to anything (such as a fluid grid), allowing for almost arbitrary geometries and straightforward global and local refinement, even during the course of a running simulation as used in publication [Pub3]. Indeed, also different discretization methodologies such as spherical harmonics are easily possible [42]. In the time domain, BIM exhibits similar freedoms, i.e.

adaptive time integration schemes can be used straightforwardly.

Furthermore, an (intrinsic⁴) interpolation from the mesh to some grid is absent. As a conse-quence, the narrow space between close objects does not need to be resolved explicitly [28].

Moreover, truly infinite domains (along one, two or all three dimensions [142,236–239] together with walls in the first two cases [240–248]) are possible. Infinite domains can be regarded as an extreme case where the fluid to particle volume ratio is zero. If it is only small (e.g. in low concentration suspensions), BIM can still have a significant performance and memory advantage over other methods (as the fluid domain is not explicitly resolved).

The boundary integral method also supports different viscosities at the inside and outside of objects without the need to keep track of internal points, which is different to e.g. Lattice-Boltzmann [212,249]. No intrinsic non-physical parameters occur by default, and the parameters often correspond directly to the ones measured in experiments. Hence, some sort of fitting or

4The smooth particle mesh Ewald method (section5.2) introduces a grid in the whole domain and a corresponding interpolation. However, this is a mere numerical trick to speed up the computation of long-range interactions rather than an intrinsic property of the method.

5 BIM: Algorithms and implementation

other conversion is not required in contrast to e.g. dissipative particle dynamics with spring-network models [33]. Finally, volume-changing objects (chapter 3) and truly solid objects (“completed double-layer boundary integral method” [132,133,168,196,197,250]) can be supported by including additional terms.

Disadvantages of BIM On the downside, an absolutely fundamental restriction of BIM is its limitation to small Reynolds numbers as it is based on the Stokes equation, meaning that any fluid inertia is necessarily absent. Furthermore, thermal fluctuations are not an intrinsic ingredient and a method to incorporate them has been only recently presented [251].

Moreover, inclusion of channel boundaries can be a bit tricky: First, direct inclusion of the walls as mesh in periodic domains leads to some flow outside of the actual channel (as the channel is embedded in the unit cell) and consequently to a (slightly) varying mean flow within the channel.

In our applications in chapter6, however, we found that this is no issue. Second, for the direct method either the velocity or the traction jump can be prescribed on the wall surface. Stipulating the velocity leads to a proper no-slip condition but an ill-conditioned Fredholm equation of the first kind, while prescribing the traction jump implies a somewhat deformable wall but good natured numerical characteristics. We usually use the latter possibility. The alternative would be to include the walls directly in the Green’s functions via explicit but complicated expressions [240–248] or the “General Geometry Ewald like method” [195,252–259].

In general, implementing boundary integral methods can be quite complicated which is at least in part due to the singular behavior of the Green’s functions, making special singularity removal procedures necessary [28,29,134,260,261]. This also applies to near-singular cases such as objects coming close to each other (which puts the above stated advantage of not having to resolve the fluid between objects somewhat into perspective).

Another undesirable property is the global aspect of the equations: Everything interacts with everything instantaneously. To prevent a quadratic scaling with the system size, advanced approaches such as the smooth particle mesh Ewald (SPME) [262] or the fast multipole method [92,263–266] must be used. The latter can also fix the troublesome parallelization property of SPME (compare section5.3).