Contact-aware simulations of particulate Stokesian suspensions
CContact-aware simulations of particulate Stokesian suspensions
Libin Lu a , Abtin Rahimian a , Denis Zorin a a Courant Institute of Mathematical Sciences, New York University, New York, NY 10003
Abstract
We present an efficient, accurate, and robust method for simulation of dense suspensions of deformableand rigid particles immersed in Stokesian fluid in two dimensions. We use a well-established boundaryintegral formulation for the problem as the foundation of our approach. This type of formulations, with ahigh-order spatial discretization and an implicit and adaptive time discretization, have been shown to beable to handle complex interactions between particles with high accuracy. Yet, for dense suspensions, verysmall time-steps or expensive implicit solves as well as a large number of discretization points are requiredto avoid non-physical contact and intersections between particles, leading to infinite forces and numericalinstability.Our method maintains the accuracy of previous methods at a significantly lower cost for dense suspen-sions. The key idea is to ensure interference-free configuration by introducing explicit contact constraintsinto the system. While such constraints are unnecessary in the formulation, in the discrete form of theproblem, they make it possible to eliminate catastrophic loss of accuracy by preventing contact explicitly.Introducing contact constraints results in a significant increase in stable time-step size for explicit time-stepping, and a reduction in the number of points adequate for stability.
1. Introduction
Particulate Stokesian suspensions of deformable and rigid particles are prevalent in nature and have manyimportant industrial applications, for example, emulsions, colloidal structures, particulate suspensions, andblood. Most of these examples are complex fluids , i.e., fluids with unusual macroscopic behavior, oftendefying a simple constitutive-law description. A major challenges in understanding the physics of complexfluids is the link between microscopic and macroscopic fluid behavior. Dynamic simulation is a powerfultool [ BK06, Naz + ] to gain insight into the underlying physical principles that govern these suspensionsand to obtain relevant constitutive relationships.Nevertheless, simulating dense suspensions of rigid and deformable particles entails many numericalchallenges, one of which is the numerical breakdown due to collision of particles, making robust treatmentof contact vital for these simulations. To this end, we present an efficient, accurate, and robust method forsimulation of dense suspensions in Stokesian fluid in D (e.g., Fig. 1), which is extendable to D . We use theboundary integral formulation to represent the flow and impose the contact-free condition as a constraint.This work focuses mainly on suspensions of rigid bodies and vesicles with high volume fractions, in whichmultiple particles are in contact or near-contact.Vesicles are closed deformable membranes suspended in a viscous medium. The dynamic deformation ofvesicles and their interaction with the Stokesian fluid play an important role in many biological phenomena.They are used to understand the properties of biomembranes [ Ghi +
11, Mis06 ] , and to simulate the motionof blood cells, in which vesicles with moderate viscosity contrast are used to model red blood cells and highviscosity contrast vesicles or rigid particles are used to model white blood cells [ Bas + ] .Boundary integral formulations offers a natural approach for accurate simulation of vesicle flows, byreducing the problem to solving equations on surfaces, and eliminating the need for discretizing changing Email addresses: [email protected] (Libin Lu), [email protected] (Abtin Rahimian), [email protected] (DenisZorin)
Preprint submitted to Elsevier October 14, 2018 a r X i v : . [ m a t h . NA ] D ec a) t = (b) t = (c) t = (d) t = Figure 1: S EDIMENTATION OF ONE HUNDRED VESICLES . Vesicles are randomly placed in a container, where only gravita-tional force is present. (a) The initial configuration. The colored tracer particles are shown for visualization purposesand are not hydrodynamically active. (b) and (c) The intermediate states; we see that as vesicles move downward theyinduce a strong back flow. (d) shows the state where the simulation was terminated. As the vesicles accumulate at thebottom of the container, many collision areas stay active. Nonetheless, the vesicles are stacked stably and without anyartifact from the collision handling.
3D volumes. However, in non-dilute suspensions, these methods are hindered by difficulties: inaccuraciesin computing near-singular integrals, and artificial force singularities caused by (non-physical) intersectionof particles. Contact situations in the Stokesian particulate flows occur frequently when the volume fractionof suspensions is high, viscosity contrast of vesicles is high, or rigid particles are present. The dynamicsof particle collision in Stokes flow are governed by the lubrication film formation and drainage, whichhas a time scale much shorter than that of the flow [ Fro + ] . Solely relying on the hydrodynamics toprevent contact requires the accurate solution of the flow in the lubrication film, which in turn entailsvery fine spatial and temporal resolution as well as expensive implicit time-stepping — imposing excessivecomputational burden as the volume fraction increases.2hile adaptive time-stepping [ QB14, QB15 ] goes a long way in maintaining stability and efficiency indilute suspensions, the time-step is determined by the closest pair of vesicles, and tends to be uniformlysmall for dense suspensions.In this work we take a different approach: we augment the governing equations with the contact con-straint. While from the point of view of the physics of the problem such a constraint is redundant, asnon-penetration is ensured by fluid forces, in numerical context it plays an important role, improving bothrobustness and accuracy of simulations. Typically, a contact law / constraint is characterized by conditions ofnon-penetration, no-adhesion as well as a mechanical complementarity condition, i.e., the contact force iszero when there is no collision. These three conditions are known as Signorini conditions in the context ofcontact mechanics or KKT conditions in the context of constrained optimization [ NW06, Wri06 ] . Contact constraints ensure that the discretized system remains intersection-free, even for relatively coarsespatial and temporal discretizations. The resulting problem is a Nonlinear Complementarity Problem (
NCP ),which we linearize and solve using an iterative method that avoids explicit construction of full matrices.We describe an implicit-explicit time-stepping scheme, adapting Spectral Deferred Correction (
SDC ) to ourconstrained setting, Section 3.2.In our approach, the minimum distance between vesicles is controlled by constraints, and is independentof the temporal resolution. While this requires solving additional auxiliary equations for the constraint forcesat every step, this additional cost is more than compensated by the ability of our method to maintain largertime-steps, and lower spatial resolutions for a given target error.For high volume fraction, our method makes it possible to increase the step size by at least an order ofmagnitude, and the simulation remain stable even for relatively coarse spatial discretizations (16 points pervesicle, versus at least 64 needed for stability without contact resolution).
We use the boundary integral formulation based on [ QB14, Rah +
10, Vee + ] ; the basic formulation usesintegral equation form of the problem and includes the effects of the viscosity contrast, fixed boundaries, aswell as deformable and rigid moving bodies. We add contact constraints to this formulation, as an inequalityconstraint on a gap function that is based on space-time intersection volume [ Har + ] . The contact force isthen parallel to the gradient of this volume with the Lagrange multiplier as its magnitude.In case of multiple simultaneous contact, this leads to a Nonlinear complementarity problem ( NCP ) forthe Lagrange multipliers, which we solve using a Newton-like matrix-free method, as a sequence of LinearComplementarity Problems (
LCP ) [ Cot +
09, Erl13 ] , solved iteratively using GMRES . The spectral Fourierbases are used for spatial discretization. For time stepping, we use semi-implicit backward Euler or semi-implicit Spectral Deferred Correction (
SDC ). Stokesian particle models are employed to theoretically and ex-perimentally investigate the properties of biological membranes [ Sac96 ] , drug-carrying capsules [ SS01 ] ,and blood cells [ NG05, Poz90 ] . There is an extensive body of work on numerical methods for Stokesianparticulate flows and an excellent review of the literature up to 2001 can be found in Pozrikidis [ Poz01 ] .Reviews of later advances can be found in [ Rah +
10, Rah +
15, Vee + ] . Here, we briefly summarize themost important numerical methods and discuss the most recent developments.Integral equation methods have been used extensively for the simulation of Stokesian particulate flowssuch as droplets and bubbles [ LH97, Loe98, RA78, ZP93 ] , vesicles [ Far +
13, Fre07, Poz90, Rah +
10, Rah + +
10, Vee +
09, ZS11, ZS13 ] , and rigid particles [ Pow93, PM87, YA75 ] . Other methods — such as phase-field approach [ Bib +
05, DZ08 ] , immersed boundary and front tracking methods [ KL10, YB12 ] , and levelset method [ Laa + ] — are used by several authors for the simulation of particulate flows.For certain flow regimes, near interaction and collision of particles has been a source of difficulty, whichwas addressed either by spatial and temporal refinement to resolve the correct dynamics (increasing the3omputational burden) or by the introduction of repulsion forces (making the time-stepping stiff).Sangani and Mo [ SM94 ] presented a framework for dynamic simulation of rigid particles with spheri-cal or cylindrical shapes, in which the lubrication forces were included directly by putting Stokes doubletsat the contact midpoint. The magnitude of lubrication force was computed using asymptotic analysis. Tomaintain the accuracy in the interaction of deformable drops, [ LH97, ZD06, Zin + ] resorted to time-steprefinement where the time step is kept proportional to particle distance d . Zinchenko and Davis [ ZD06 ] ,Zinchenko et al. [ Zin + ] keep the time step proportional to (cid:112) d . Loewenberg and Hinch [ LH97 ] adjustboth the grid spacing around the contact region and the time-step to be proportional to d . Freund [ Fre07 ] resorted to repulsion force to avoid contact in a D particulate flow. In a later work for D , Freund and coau-thors [ Zha + ] observed that significantly larger repulsion force density are needed in three dimensions,as the total repulsion force is distributed over a smaller region, when measured as a fraction of the totalsurface area / length. Consequently, they used a purely kinematic collision handing, in which, after eachtime-step, the intersecting points are moved outside.Quaife and Biros [ QB15 ] applies adaptive time-stepping and backtracking to resolve collisions. Similarly,Ojala and Tornberg [ OT14 ] present an interesting integral equation method for the flow of droplets in twodimensions with a specialized quadrature scheme for accurate near-singular evaluation enabling simulationof flows with close to touching particles. While methods using adaptivity both in space and time are themost robust and accurate, they incur excessive cost as means of collision handling. Related work on contact response.
A broad range of methods were developed for collision detection and re-sponse. While the work in contact mechanics often focuses on capturing the physics of the contact correctly(e.g., taking into account friction effects), the work in computer graphics literature emphasizes robustnessand efficiency. In our context, robustness and efficiency are particularly important, as we aim to modelvesicle flows with high density and large numbers of particles. Physical correctness has a somewhat differ-ent meaning: as we know that if the forces and surfaces in the system are resolved with high accuracy, thecontacts would not occur, our primary emphasis is on reducing the impact of the artificial forces associatedwith contacts on the system.There is an extensive literature on contact handling in computational contact mechanics mainly in thecontext of
FEM mechanical and thermal analysis [ FW05, JJ87, KW02, Pus +
08, Tur +
09, Wri95, Wri06 ] .Wriggers [ Wri95, Wri06 ] presents in-depth reviews of the contact mechanics framework. The works incontact mechanics literature can be categorized based on their ability in handling large deformations and / ortangential friction. In some of the methods, to simplify the problem, small deformation assumption is used topredefine the active part of the boundary as well as to align the FEM mesh. Numerical methods for contactresponse can be categorized as (i) penalty forces, (ii) impulse / kinematic responses, and (iii) constraintsolvers.From algorithmic viewpoint, contact mechanics methods in FEM include: (i) Node-to-node methodswhere the contact between nodes is only considered. The
FEM nodes of contacting bodies need to alignedand therefore this method is only applicable to small deformation. (ii) Node-to-surface methods check thecollision between predefined set of nodes and segments. Similar to node-to-node methods, these meth-ods can only handle small deformations. (iii) Surface-to-surface methods, where the contact constraint isimposed in weak form. In contrast to the two previous class of algorithms, methods in this class are capa-ble of handling large deformations. Mortar Method is well-known within this class of algorithms [ FW05,KW02, Pus +
08, Tur + ] . The Mortar Method was initially developed for connecting different non-matchingmeshes in the domain decomposition approaches for parallel computing, e.g., [ Pus04 ] .In these methods, no-penetration is either enforced as a constraint using a Lagrangian (identified withthe contact pressure) or penalty force based on a gap function. To the best of our knowledge, for contactmechanics problems, a signed distance between geometric primitives is used as the gap function, in contrastto our approach where we use space-time interference volume.Fischer and Wriggers [ FW05 ] present a frictionless contact resolution framework for D finite defor-mation using Mortar Method using penalty force or Lagrange multiplier. Tur et al. [ Tur + ] use similarmethod for frictional contact in D . Puso et al. [ Pus + ] use Mortar Method for large deformation contact4sing quadratic element.Our problem has similarities to large-deformation frictionless contact problems in contact mechanics.An important difference however, is the presence of fluid, which plays a major role in contact response.Application of boundary integral methods in contact mechanics is rather limited compared to the FEM methods [ Eck +
99, Gun04 ] . Eck et al. [ Eck + ] used Boundary Element Method for the static contactproblem where Coulomb friction is presented. Gun [ Gun04 ] solved static problem with load increment andcontact constraint on displacement and traction.In computer graphics literature, a set of commonly used and efficient methods are based on [ Pro97 ] ,a method for the collision handling of mass-spring cloth models. To ensure that the system remainsintersection-free, zones of impact are introduced and rigid body motion is enforced in each zone of impact;while this method works well in practice, its effects on the physics of the objects are difficult to quantify.Penalty methods are common due to the ease in their implementation, but suffer from time-steppingstiffness and / or the lack of robustness. Baraff and Witkin [ BW98 ] uses implicit time-stepping coupled withrepulsion force equal to the variation of the quadratic constraint energy with respect to control vertices.Soft collisions are handled by the introduction of damped spring and rigid collisions are enforced by mod-ification to the mass matrix. Faure et al. [ Fau + ] introduced Layered Depth Images to allow efficientcomputation of the collision volumes and their gradients using GPU s. A penalty force proportional to thegradient is used to resolve collisions. However, the stiffness of the repulsion force varies greatly (from 10 to 10 ) in their experiments. To address these difficulties, Harmon et al. [ Har + ] present a frameworkfor robust simulation of contact mechanics using penalty forces through asynchronous time-stepping, albeitat a significant computational cost. Alternatively, one can view collision response as an instantaneous reac-tion (an impulse), i.e., an instantaneous adjustment of the velocities. However, such adjustments are oftenproblematic in the case of multiple contacts, as these may lead to a cyclic “trembling” behavior.Our method belongs to a large family of constraint-based methods, which are increasingly the standardapproach to contact handling in graphics. This set of methods meets our goals of providing robustnessand improving efficiency of contact response, while minimizing the impact on the physics of the system (asvirtual forces associated with constraints do not do work).Duriez et al. [ Dur + ] start from Signorini’s law and derive the contact force formulation. The resultingequation is an LCP that is solved by Gauss–Seidel like iterations, sequentially resolving contacts until reach-ing the contact free state. Harmon et al. [ Har + ] focus on robust treatment of collision without simulationartifacts. To enforce the no-collision constraint, this work uses an impulse response that gives rise to an LCP problem for its magnitude. To reduce the computational cost, the
LCP solution (the Lagrange multiplier) isapproximated by solving a linear system. Otaduy et al. [ Ota + ] uses a linear approximation to contactconstraints and a semi-implicit discretization, solving a mixed LCP problem at each iteration.Our approach is directly based on [ Har + ] and is closest to [ All + ] , in which the intersection volumeand its gradient with respect to control vertices are computed at the candidate step. The non-collisionis enforced as a constraint on this volume, which lead to a much smaller system compared to distanceformulation between geometric primitives. The constrained formulation leads to an LCP problem. Harmonet al. [ Har + ] assumes linear trajectory between edits and define space-time interference volume and usesit as a gap function and we use similar formulation to define the interference volume. In Table 1 we list symbols and operators frequently used in this paper. Throughout this paper, lower caseletters refer to scalars, and lowercase bold letters refer to vectors. Discretized quantities are denoted bysans serif letters.
2. Formulation
We consider the Stokes flow with N v vesicles and N p rigid particles suspended in a Newtonian fluid whichis either confined or fills the free space, Fig. 2. In Stokesian flow, due to high viscosity and / or small lengthscale, the ratio of inertial and viscous forces (The Reynolds number) is small and the fluid flow can be5 ymbol Definition γ i The boundary of the i th vesicle γ ∪ i γ i ν i The viscosity contrast µ Viscosity of the ambient fluid µ i Viscosity of the fluid inside i th vesicle σ Tension χ Shear rate (cid:36) i The domain enclosed by γ i (cid:36) ∪ i (cid:36) i d Separation distance of particles d m Minimum separation distance f σ Tensile force f b Bending force f c Collision force h Distance between two discretizationpoints (on each surface) Symbol Definition J Jacobian of contact volumes V n Unit outward normal u Velocity u ∞ The background velocity field V Contact volumes X Coordinate of a (Lagrangian) point on asurface (cid:71)
Stokes Single-layer operator (cid:84)
Stokes Double-layer operator
LCP
Linear Complementarity problem
NCP
Nonlinear Complementarity Problem
SDC
Spectral Deferred Correction
STIV
Space-Time Interference Volumes LI Locally-implicit time-stepping
CLI
Locally-implicit constrained time-stepping GI Globally-implicit time-stepping
Table 1: I NDEX OF FREQUENTLY USED SYMBOLS , OPERATORS , AND ABBREVIATIONS . described by the incompressible Stokes equation − µ ∆ u ( x ) + ∇ p ( x ) = F ( x ) , and ∇· u ( x ) = ( x ∈ Ω) , (2.1) F ( x ) = (cid:90) γ f ( X ) δ ( x − X ) d s ( X ) , (2.2)where f is the surface density of the force exerted by the vesicle’s membrane on the fluid. Ω denotesthe fluid domain of interest with Γ as its enclosing boundary (if present) and µ denoting the viscosityof ambient fluid. If Ω is multiply-connected, its interior boundary consists of K smooth curves denotedby Γ , . . . , Γ K . The outer boundary Γ encloses all the other connected components of the domain. Theboundary of the domain is then denoted Γ : = (cid:83) k Γ k . We use x to denote an Eulerian point in the fluid( x ∈ Ω ) and X a Lagrangian point on the vesicles or rigid particles. We let γ i denote the boundary of the i th vesicle ( i =
1, . . . , N v ), (cid:36) i denote the domain enclosed by γ i , µ i denote viscosity of the fluid inside thatvesicle, and γ : = (cid:83) i γ i . Equation (2.1) is valid for x ∈ (cid:36) i by replacing µ with µ i .There are rigid particles suspended in the fluid domain. We denote the boundary of the j th rigid particleby π j ( j =
1, . . . , N p ) and let π : = (cid:83) j π j . The governing equations are augmented with the no-slip boundarycondition on the surface of vesicles and particles u ( X , t ) = X t ( X ∈ γ ∪ π ) , (2.3)where X t : = ∂ X ∂ t is the material velocity of point X on the surface of vesicles or particles. The velocity on thefixed boundaries is imposed as a Dirichlet boundary condition u ( x ) = U ( x ) ( x ∈ Γ) . (2.4)We assume that the vesicle membrane is inextensible, i.e., X s · u s = ( X ∈ γ ) , (2.5)where the subscript “ s ” denotes differentiation with respect to the arclength on the surface of vesicles.Rigid particles are typically force- and torque-free. However, surface forces may be exerted on them dueto a constraint, e.g., the contact force f c , which we will define later. In this case, the force F π j and torque6 Γ Γ π π γ γ Collision n n
Figure 2: S CHEMATIC . The flow domain Ω (gray shaded area) with boundary Γ k ( k =
0, . . . , K ) . Vesicles and rigidparticles are suspended in the fluid. The vesicle boundaries are denoted by γ i ( i =
1, . . . , N v ) and the rigid particles(checkered pattern) are denoted by π j ( j =
1, . . . , N p ) . The outward normal vector to the boundaries is denotes by n . The dotted lines around boundaries denote the prescribed minimum separation distance for each of them. Theminimum separation distance is a parameter and can be set to zero. In this schematic, vesicle γ and particle π aswell as vesicle γ and boundary Γ are in contact. The slices of the space-time intersection volumes at the currentinstance are marked by . L π j exerted on the j th particle are the sum of such terms induced by constraints F π j = , or F π j = (cid:90) π j f c ( X ) d s ( j =
1, . . . , N p ) , (2.6) L π j = , or L π j = (cid:90) π j ( X − c π j ) · f c ⊥ ( X ) d s ( j =
1, . . . , N p ) , (2.7)where c π j is the center of mass for π j and for vector f = ( f , f ) , f ⊥ : = ( f , − f ) . It is known [ Fro +
13, Nem + ] that the exact solution of equations of motion, Eqs. (2.1), (2.3), and (2.4),keeps particles apart in finite time due to formation of lubrication film. Thus, it is theoretically sufficient tosolve the equations with an adequate degree of accuracy to avoid any problems related to overlaps betweenparticles. Nonetheless, achieving this accuracy for many types of flows (most notably, flows with highvolume fraction of particles or with complex boundaries) may be prohibitively expensive.With inadequate computational accuracy particles may collide with each other or with the boundariesand depending on the numerical method used, the consequences of this may vary. For methods basedon integral equations the consequences are particularly dramatic, as overlapping boundaries may lead todivergent integrals. To address this issue, we augment the governing equations with a contact constraint,formally written as V ( S , t ) ≤
0, (2.8)7here S = Γ ∪ γ ∪ π denotes the boundary of the domain and all particles. The function V is chosen in away that V > S are at a distance less than a user-specified constant d m .Function V may be a vector-valued function, for which the inequality is understood component-wise. Thisconstraint ensures that the suspension remains contact-free independent of numerical resolution.For the constraint function V , in addition to the basic condition above, we choose a function that satisfiesseveral additional criteria:(i) it introduces a relatively small number of additional constraints, and(ii) when the function is discretized, no contacts are missed even for large time steps.To clarify the second condition, suppose we have a small particle rapidly moving towards a planar boundary.For a large time step, it may move to the other side of the boundary in a single step, so any condition thatconsiders an instantaneous quantity depending on only the current position is likely to miss the contact.To this end, we use the Space-Time Interference Volumes ( STIV ) from Harmon et al. [ Har + ] to definethe function V as the area in space-time swept by the intersecting segments of boundary over time. To bemore precise, for each point X on the boundary S , consider a trajectory X ( τ ) , between a time t , for whichthere are no collisions, and a time t . Points X ( τ ) define a deformed boundary S ( τ ) for each τ . For eachpoint X , we define τ I ( X ) , t ≤ τ I ≤ t , to be the first instance for which this point comes into contact witha different point of S ( τ I ) . We let I ( X ) denote the point on S ( τ I ) that comes into contact with X ( τ I ) , i.e., X ( τ I ) = I ( X )( τ I ) . For points which never come into contact with other points, we set τ I = t . Then, thespace-time interference volume for D problems is defined as V ( S , t ) = (cid:90) S ( t ) (cid:90) t τ I ( X ) (cid:104) + (cid:0) X t ( τ ) · n ( X , τ ) (cid:1) (cid:105) / d τ d s (2.9)where n ( X , τ ) denotes the normal to S ( τ ) at X ( τ ) . The integration is over the initial boundary S ( t ) , andwe use the fact that the surface is inextensible and the surface metric does not change. Note that, the termvolume is a misnomer because Eq. (2.9) is the surface area swept by the intersecting segments of boundary.An important property of this choice of function, compared to, e.g., a space intersection volume, is that foreven a very thin object moving at high velocity, it will be proportional to the time interval t − t .We consider each connected component of this volume as a separate volume, and impose an inequalityconstraint on each; while keeping a single volume is in principle equivalent, using multiple volumes avoidcertain undesirable effects in discretization [ Har + ] . Thus, V ( S , t ) is a vector function of time-dependentdimension, with one component per active contact region.Depending on the context, we may omit the dependence of V on S and write V ( t ) as the contact volumefunction or V ( γ i , t ) for the sub-vector of V ( S , t ) involving surface γ i .In practice, it is desirable to control the minimal distance between particles. Therefore, we define aminimum separation distance d m ≥ d m distance from each other; as shown in Fig. 2. The contact volume with minimumseparation distance is calculated with the surface displaced by d m , i.e., the time t I is obtained not from thefirst contact with S ( τ ) but rather the displaced surface S ( τ ) + d m n ( τ ) . Maintaining minimum separationdistance — rather than considering pure contact only — eliminates of potentially expensive computation ofnearly singular integrals close to the surface and improves the accuracy in semi-explicit time-stepping. We use the Lagrange multiplier method (e.g., [ Wri06 ] ) to add contact constraints to the system. While itis computationally more expensive than adding a penalty force for the constraint (effectively, an artificialrepulsion force), it has the advantage of eliminating the need of tuning the parameters of the penalty forceto ensure that the constraint is satisfied and keeping nonphysical forces introduced into the system to the8inimum required for maintaining the desired separation. The constrained system can be written asmin (cid:90) Ω (cid:18) µ ∇ u · ∇ u − u · F (cid:19) d A , (2.10)subject to: ∇· u ( x ) = ( x ∈ Ω) , X s · u s = ( X ∈ γ ) , V ( S , t ) ≥ (cid:76) ( u , p , σ , λ ) = (cid:90) Ω (cid:18) µ ∇ u · ∇ u − u · F − p ∇· u (cid:19) d A + (cid:90) γ σ X s · u s d s + λ · V . (2.11)The first-order optimality ( KKT ) conditions yield the following modified Stokes equation, along with theconstraints listed in Eq. (2.10): − µ ∆ u + ∇ p = F (cid:48) , (2.12) F (cid:48) = F + (cid:90) γ f σ δ ( x − X ) d s + (cid:90) S f c δ ( x − X ) d s , (2.13) f σ = − ( σ X s ) s , (2.14) f c = ∇ V T λ , (2.15) λ ≥
0, (2.16) λ · V =
0, (2.17)where the last condition is the complementarity condition — either an equality constraint is active ( V i = ) or its Lagrange multiplier is zero. As we will see in the next section and based on Eq. (2.13), the collisionforce f c is added to the traction jump across the vesicle’s interface. For rigid particles, the contact forceinduces force and torque on each particle — as given in Eqs. (2.6) and (2.7).It is customary to combine V ≥ λ ≥
0, and λ · V =
0, into one expression and write0 ≤ V ( t ) ⊥ λ ≥
0, (2.18)where " ⊥ " denotes the complementarity condition. These ensure that the Signorini conditions introducedin Section 1 are respected: contacts do not produce attraction force ( λ ≥
0) and the constraint is active ( λ nonzero) if and only if V ( t ) is zero. Furthermore, it follows from V ( t ) · λ = f c isperpendicular to the velocity and therefore it respects the principal of virtual work and does not add to orremove from the system’s energy. Following the standard approach of potential theory [ Pow93, Poz92 ] , one can express the solution of theStokes boundary value problem, Eq. (2.12), as a system of singular integro-differential equations on allimmersed and bounding surfaces.The Stokeslet tensor G , the Stresslet tensor T , and the Rotlet R are the fundamental solutions of theStokes equation given by G ( r ) = πµ (cid:130) − log (cid:107) r (cid:107) I + r ⊗ r (cid:107) r (cid:107) (cid:140) , (2.19) T ( r ) = πµ r ⊗ r ⊗ r (cid:107) r (cid:107) , (2.20) R ( r ) = πµ r ⊥ (cid:107) r (cid:107) . (2.21)9he solution of Eq. (2.12) can be expressed by the combination of single- and double-layer integrals.We denote the single-layer integral on the vesicle surface γ i by (cid:71) γ i [ f ]( x ) : = (cid:90) γ i G ( x − Y ) · f ( Y ) d s , (2.22)where f is an appropriately defined density. The double-layer integral on a surface S (a vesicle, a rigidparticle, or a fixed boundary) is (cid:84) S [ q ]( x ) : = (cid:90) S n ( Y ) · T ( x − Y ) · q ( Y ) d s , (2.23)where n denotes the outward normal (as shown in Fig. 2) to the surface S , and q is appropriately defineddensity. When the evaluation point x is on the integration surface, Eq. (2.22) is a singular integral, andEq. (2.23) is interpreted in the principal value sense.Due to the linearity of the Stokes equations, as formulated in [ QB14, Rah + ] , the velocity at a point x ∈ Ω can be expressed as the superposition of velocities due to vesicles, rigid particles, and fixed boundaries α u ( x ) = u ∞ ( x ) + u γ ( x ) + u π ( x ) + u Γ ( x ) , x ∈ Ω , α = x ∈ Ω \ γ , ν i x ∈ (cid:36) i , ( + ν i ) / x ∈ γ i , (2.24)where u ∞ ( x ) represent the background velocity field (for unbounded flows) and ν i = µ i /µ denotes theviscosity contrast of the i th vesicle. The velocity contributions from vesicles, rigid particles, and fixed bound-aries each can be further decomposed into the contribution of individual components u γ ( x ) = N v (cid:88) i = u γ i ( x ) , u π ( x ) = N p (cid:88) j = u π j ( x ) , u Γ ( x ) = K (cid:88) k = u Γ k ( x ) . (2.25)To simplify the representation, we introduce the complementary velocity for each component. For the i th vesicle, it is defined as ¯ u γ i = α u − u γ i . The complementary velocity is defined in a similar fashion for rigidparticles as well as components of the fixed boundary. The velocity induced by the i th vesicle is expressed as an integral [ Poz92 ] : u γ i ( x ) = (cid:71) γ i [ f ]( x ) + ( − ν i ) (cid:84) γ i [ u ]( x ) ( x ∈ Ω) , (2.26)where the double-layer density u is the total velocity and f is the traction jump across the vesicle membrane.Based on Eq. (2.13), the traction jump is equal to the sum of bending, tensile, and collision forces f ( X ) = f b + f σ + f c = − κ b X ssss − ( σ X s ) s + ∇ V T λ ( X ∈ γ ) , (2.27)where κ b is the membrane’s bending modulus. The tensile force f σ = ( σ X s ) s is determined by the localinextensibility constraint, Eq. (2.5), and tension σ is its Lagrangian multiplier.Note that Eq. (2.26) is the contribution from each vesicle to the velocity field. To obtain an equation forthe interfacial velocity, Eq. (2.26) is to be substituted into Eq. (2.24): ( + ν i ) u ( X ) = ¯ u γ i ( X ) + (cid:71) γ i [ f ]( X ) + ( − ν i ) (cid:84) γ i [ u ]( X ) ( X ∈ γ i ) , (2.28)subject to local inextensibility X s · u s = ( X ∈ γ i ) . (2.29)10 .3.2. The contribution from the fixed boundaries. The velocity contribution from the fixed boundary can beexpressed as a double-layer integral [ Pow93 ] along Γ . The contribution of the outer boundary Γ is u Γ ( x ) = (cid:84) Γ [ η ]( x ) ( x ∈ Ω) , (2.30)where η is the density to be determined based on boundary conditions. Substituting Eq. (2.30) intoEq. (2.24) and taking its limit to a point on Γ and using the Dirichlet boundary condition, Eq. (2.4), weobtain a Fredholm integral equations for the density η U ( x ) − ¯ u Γ ( x ) = − η ( x ) + (cid:84) Γ [ η ]( x ) ( x ∈ Γ ) .However, this equation is rank deficient. To render it invertible, the equation is modified following [ Poz92 ] : U ( x ) − ¯ u Γ ( x ) = − η ( x ) + (cid:84) Γ [ η ]( x ) + (cid:78) Γ [ η ]( x ) ( x ∈ Γ ) , (2.31)where the operator (cid:78) Γ is defined as (cid:78) Γ [ η ]( x ) = (cid:90) Γ [ n ( x ) ⊗ n ( y )] · η ( y ) d s . (2.32)For the enclosed boundary components Γ k ( k > ) to eliminate the double-layer nullspace we need toinclude additional Stokeslet and Rotlet terms u Γ k ( x ) = (cid:84) Γ k [ η k ]( x ) + G ( x − c Γ k ) · F Γ k + R ( x − c Γ k ) L Γ k , ( k =
1, . . . , K ; x ∈ Ω) , (2.33)where c Γ k is a point enclosed by Γ k , F Γ k is the force exerted on Γ k , and L Γ k is the torque: F Γ k = | Γ k | (cid:90) Γ k η k d s , L Γ k = | Γ k | (cid:90) Γ k ( X − c Γ k ) · η ⊥ k d s , (2.34)where | Γ k | denotes the perimeter of Γ k . Taking the limit at points of the surface, leads to the followingintegral equation: U ( x ) − ¯ u Γ k ( x ) = − η k ( x ) + (cid:84) Γ k [ η k ]( x ) + G ( x − c Γ k ) · F Γ k + R ( x − c Γ k ) L Γ k ( x ∈ Γ k ) . (2.35)Equations (2.34) and (2.35) are a complete system for double-layer densities η k , forces F Γ k , and torques L Γ k on each enclosed surface Γ k . The formulation for rigid particles is very similar to that of fixedboundaries, except the force and torque are known (cf. Eqs. (2.6) and (2.7)). The velocity contributionfrom the j th rigid particle is u π j ( x ) = (cid:84) π j [ ζ j ]( x ) + G ( x − c π j ) · F π j + R ( x − c π j ) L π j , (2.36)Where F π j , L π j are, respectively, the known net force and torque exerted on the particle and ζ j is the unknowndensity.Let U π j and ω π j be the translational and angular velocities of the j th particle; then we obtain the followingintegral equation for the density ζ j from the limit of (2.36): U π j + ω π j ( X − c π j ) ⊥ − ¯ u π j ( X ) = − ζ j ( X ) + (cid:84) π j [ ζ j ]( X ) + G ( X − c π j ) · F π j + R ( X − c π j ) L π j . (2.37)where F π j = | π j | (cid:90) π j ζ j d s , L π j = | π j | (cid:90) π j ( Y − c π j ) · ζ ⊥ j d s (2.38)where c π j is the center of j th rigid particle. Equations (2.37) and (2.38) are used to solve for the unknowndensities ζ j as well as the unknown translational and angular velocities of each particle. Note that theobjective of Eqs. (2.34) and (2.38) is remove the null space of the double-layer and therefore the left-hand-side (i.e. the image of the null function) can be chosen rather arbitrarily.11 .4. Formulation summary The formulae outlined above govern the evolution of the suspension. The flow constituents are hydrody-namically coupled through the complementary velocity (i.e., the velocity from all other constituents). Giventhe configuration of the suspension, the unknowns are: • Velocity u ( X ) and tension σ of vesicles’ interface determined by Eqs. (2.27–2.29). The velocity isintegrated for the vesicles’ trajectory using Eq. (2.3). • The double-layer density on the enclosing boundary η as well as the double-layer density η k ( k =
1, . . . , K ) , force F Γ k , and torque L Γ k on the interior boundaries determined by Eqs. (2.31), (2.34),and (2.35). Note that the collision constraint does not enter the formulation for the fixed boundariesand when a particle collides with a fixed boundary, the collision force is only applied to the particle.The unknown force and torque above can be interpreted as the required force to keep the boundaryin place. • Translational U π j and angular ω π j velocities of rigid particles ( j =
1, . . . , N p ) as well as double-layerdensities ζ j on their boundary determined by Eqs. (2.37) and (2.38). Where the force and torque areeither zero or determined by the collision constraint Eqs. (2.6) and (2.7).This system is constrained by Signorini ( KKT ) conditions for the contact, Eq. (2.18), which is used tocompute λ , the strength of the contact force.In the referenced equations above, the complementary velocity is combination of velocities given inEqs. (2.26), (2.30), (2.33), and (2.36).
3. Discretization and Numerical Methods
In this section, we describe the numerical algorithms required for solving the dynamics of a particulateStokesian suspension. We use the spatial representation and integral schemes in [ Rah + ] . We also adaptthe spectral deferred correction time-stepping in [ QB15, QB16 ] to the local implicit time-stepping schemes.Furthermore, we use piecewise-linear discretization of curves to calculate the space-time contact volume V ( γ , t ) , Eq. (2.9), as in [ Har + ] . To solve the complementarity problem resulting from the contact con-straint, we use the minimum-map Newton method discussed in [ Erl13 ] .The key difference, compared to previous works is that at every time step instead of solving a linearsystem we solve a nonlinear complementarity problem ( NCP ). The
NCP s are solved iteratively, using a Lin-ear Complementarity Problem (
LCP ) solver. We refer to these iterations as contact-resolving iterations, incontrast to the outer time-stepping iterations.For simplicity, we describe the scheme for a system including vesicles only, without boundaries or rigidparticles. Adding these requires straightforward modifications to the equations. In the following sections,we will first summarize the spatial discretization, then discuss the
LCP solver, and close with the timediscretization with contact constraint.
All interfaces are discretized with N uniformly-spaced discretization points [ Rah + ] . The number of pointson each curve is typically different but for the sake of clarity we use N . The distance between discretizationpoints does not change with time over the curves due to rigidity of particles or the local inextensibilityconstraint for vesicles. Let X ( s ) , with s ∈ ( L ] , be a parametrization of the interface γ i ( or π j ) , and let { s k = kL / N } Nk = be N equally spaced points in arclength parameter, and X k : = X ( s k ) the correspondingmaterial points. High-order discretization for force computation.
We use the Fourier basis to interpolate the positions andforces associated with sample points, and
FFT to calculate the derivatives of all orders to spectral accuracy.12e use the hybrid Gauss-trapezoidal quadrature rules of [ Alp99 ] to integrate the single-layer potential for X ∈ γ i (cid:71) γ i [ f ]( X ) ≈ G γ i [ f ]( X ) : = N + M (cid:88) (cid:96) = w (cid:96) G ( X − Y (cid:96) ) · f ( Y (cid:96) ) , (3.1)where w (cid:96) are the quadrature weights given in [ Alp99, Table 8 ] and Y (cid:96) are quadrature points. For X ∈ γ i ,the linear operator G γ i is a matrix, we denote the single-layer potential on γ i as G γ i f .The double-layer kernel n ( Y ) · T ( X − Y ) in Eq. (2.23) is non-singular in two dimensionslim γ i (cid:51) X → Y n ( Y ) · T ( X − Y ) = − κ πµ t ⊗ t ,where t denotes the tangent vector. Therefore, a simple uniform-weight composite trapezoidal quadraturerule has spectral accuracy in this case. Similar to the single-layer case, we denote the discrete double-layerpotential on γ i by T γ i u . We use the nearly-singular integration scheme described in [ QB14 ] to maintainhigh integration accuracy for particles closely approaching each other. Piecewise-linear discretization for constraints..
While the spectral spatial discretization is used for mostcomputations, it poses a problem for the minimal-separation constraint discretization. Computing para-metric curve intersections, an essential step in the
STIV computation, is relatively expensive and difficult toimplement robustly, as this requires solving nonlinear equations. We observe that the impact of the sep-aration distance on the overall accuracy is low in most situations, as explored in Section 4. Thus, ratherthan enforcing the constraint as precisely as allowed by the spectral discretization, we opt for a low-order,piecewise-linear discretization in this case, and use an algorithm that ensures that at least the target minimalseparation is maintained, but may enforce a higher separation distance.For the purpose of computing
STIV and its gradient, we use L ( X , r ) , the piecewise-linear interpolant of r times refinement of points — the upsampled points correspond to arclength values with spacing L / ( N r ) ,with r determined dynamically.For discretized computations, we set the separation distance to ( + α ) d m , where d m is the targetminimum separation distance. We choose r such that (cid:107) L ( X , r ) − X (cid:107) ∞ < α d m . Our NCP solver, describedbelow, ensures that the separation between L ( X , r ) is ( + α ) d m at the end of a single time step. We choose α = r = α require more refinement, but enforcethe constraint more accurately.At the end of the time step, the minimal-separation constraint ensures that L ( X , r )( s ) , for any s , is atleast at the distance ( + α ) d m from a possible intersection if its trajectory is extrapolated linearly. Bycomputing the upper bounds on the difference between the X ( s ) and L ( X , r )( s ) at the beginning of thetime step, and interpolated velocities, we obtain a lower bound on the actual separation distance d (cid:48) forthe spectral surface X ( s ) . If d (cid:48) < d m , we increase r , and repeat the time step. As the piecewise linearapproximation converges to the spectral boundary X , and so do the interpolated velocities. In practice, wehave not observed a need for refinement for our choice of α .For the piecewise-linear discretization of curves, the space-time contact volume V ( γ , t ) , Eq. (2.9), and itsgradient are calculated using the definitions and algorithms in [ Har + ] . Given a contact-free configurationand a candidate configuration for the next time step, we calculate the discretized space-time contact volumeas the sum of edge-vertex contact volumes V = (cid:80) k V k ( e , X ) . We use a regular grid of size proportional tothe average boundary spacing to quickly find potential collisions. For all vertices and edges, the boundingbox enclosing their initial (collision-free) and final (candidate position) location is formed and all the gridboxes intersecting that box are marked. When the minimal separation distance d m >
0, the bounding boxis enlarged by d m . For each edge-vertex pair e ( X i , X i + ) and X j , we solve a quartic equation to find theirearliest contact time τ I assuming linear trajectory between initial and candidate positions. We calculate theedge-vertex contact volume using Eq. (2.9): V k ( e , X ) = ( t − τ I )( + ( X ( τ I ) · n ( τ I )) ) / | e | , (3.2)13here n ( τ I ) is the normal to the edge e ( τ I ) . For each edge-vertex contact volume, we calculate the gradientrespect to the vertices X i , X i + and X j , summing over all the edge-vertex contact pairs we get the totalspace-time contact volume and gradient. Our temporal discretization is based on the locally-implicit time-stepping scheme in [ Rah + ] — adaptingthe Implicit–Explicit ( IMEX ) scheme [ Asc + ] for interfacial flow — in which we treat intra-particle in-teractions implicitly and inter-particle interactions explicitly. We combine this method with the minimal-separation constraint. We refer to this scheme as constrained locally-implicit ( CLI ) scheme. For comparisonpurposes, we also consider the same scheme without constraints ( LI ) and the globally semi-implicit ( GI )scheme, where all interactions treated implicitly [ Rah + ] . From the perspective of boundary integral for-mulation, the distinguishing factor between LI and CLI is the extra traction jump term due to collision.Schemes LI / CLI and GI differ in their explicit or implicit treatment of the complementary velocities.While treating the inter-vesicle interactions explicitly may result in more frequent violations of minimal-separation constraint, we demonstrate that in essentially all cases the CLI scheme is significantly moreefficient than both the GI and LI schemes because these schemes are costlier and require higher spatial andtemporal resolution to prevent collisions.We consider two versions of the CLI scheme, a simple first-order Euler scheme and a spectral deferredcorrection version. A first-order backward Euler
CLI time stepping formulation for Eq. (2.28) is1 + ν i u + i = ¯ u γ i + G γ i f i ( X + i , σ + i , λ + ) + ( − ν i ) T γ i u + i , (3.3) X i , s · u + i , s =
0, (3.4) f i ( X + i , σ + i , λ + ) = − κ b X + i , ssss − ( σ + X i , s ) s + (( ∇ V + ) T λ + ) i , (3.5)0 ≤ V ( γ ; t + ) ⊥ λ + ≥
0, (3.6)where the implicit unknowns to be solved for at the current step are marked with superscript “ + ”. Theposition and velocity of the points of i th vesicle are denoted by X i , u + i = ( X + i − X i ) / ∆ t , and f i is the tractionjump on the i th vesicle boundary. V ( γ ; t + ) is the STIV function.
We use spectral deferred correction (
SDC ) method [ QB15, QB16 ] to geta better stability behavior compared to the basic backward Euler scheme described above. We use SDC bothfor LI and CLI time-stepping. To obtain the
SDC time-stepping equations, we reformulate Eq. (2.3) as aPicard integral X ( t n + ) = X ( t n ) + (cid:90) t n + t n u ( X ( τ ) , σ ( τ ) , λ ( τ )) d τ , (3.7)where the velocity satisfies Eqs. (2.28) and (3.3). In the SDC method, the integral in Eq. (3.7) is firstdiscretized with p + p provisional positions (cid:101) X corresponding to times τ i in the interval [ t n , t n + ] ; t n = τ < · · · < τ p = t n + . Provisional tensions (cid:101) σ andprovisional (cid:101) λ are defined similarly. The SDC method iteratively corrects the provisional positions (cid:101) X with theerror term (cid:101) e , which is solved using the residual (cid:101) r resulting from the provisional solution as defined below.The residual is given by: (cid:101) r ( τ ) = (cid:101) X ( t n ) − (cid:101) X ( τ ) + (cid:90) τ t n (cid:101) u ( θ ) d θ . (3.8)After discretization, we use (cid:101) X w , m , to denote the provisional position at m th Gauss-Lobatto point after w SDC passes. The error term (cid:101) e w , m denotes the computed correction to obtain m th provisional position in w th SDC correction iteration is defined by (cid:101) X w , m = (cid:101) X w − m + (cid:101) e w , m , (cid:101) σ w , m = (cid:101) σ w − m + (cid:101) e w , m σ , (3.9) (cid:101) λ w , m = (cid:101) λ w − m + (cid:101) e w , m λ .Setting (cid:101) X m to zero, the first SDC pass is just backward Euler time stepping to obtain nontrivial provisionalsolutions. Beginning from the second pass, we solve for the error term as corrections.Denote ( α I − ( − ν ) T (cid:101) γ w − m ) by D (cid:101) γ w − m . Following [ QB15, QB16 ] , we solve the following equation forthe error term: α (cid:101) e w , m − (cid:101) e w , m − ∆ τ = D (cid:101) γ w − m (cid:130)(cid:101) r w − m − (cid:101) r w − m − ∆ τ (cid:140) + G (cid:101) γ w − m f ( (cid:101) e w , m , (cid:101) e w , m σ , (cid:101) e w , m λ )+( − ν ) T (cid:101) γ w − m (cid:130) (cid:101) e w , m − (cid:101) e w , m − ∆ τ (cid:140) .(3.10)Eq. (3.10) is the identical to Eq. (3.3), except the right-hand-side for Eq. (3.10) is obtained from the residualwhile the right-hand-side for Eq. (3.3) is the complementary velocity. The residual (cid:101) r w , m is obtained using adiscretization of Eq. (3.8): (cid:101) r w , m = (cid:101) X w ,0 − (cid:101) X w , m + p (cid:88) l = w l , m (cid:101) u w , l . (3.11)where w l , m are the quadrature weights for Gauss-Lobatto points, whose quadrature error is (cid:79) (∆ t p − ) . Inaddition to the SDC iteration, Eq. (3.10), we also enforce the inextensibility constraint (cid:101) X w − ms · (cid:101) u w , ms =
0, (3.12)and the contact complementarity 0 ≤ V ( (cid:101) γ w , m ) ⊥ (cid:101) e w , m λ ≥
0. (3.13)In evaluating the residuals using Eq. (3.11), provisional velocities are required. In the GI scheme [ QB16 ] ,all the interactions are treated implicitly and given provisional position (cid:101) X w , m , the provisional velocities areobtained by evaluating (cid:101) u w , m = D − (cid:101) γ w , m (cid:16) G (cid:101) γ w , m f ( (cid:101) X w , m , (cid:101) σ w , m , (cid:101) λ w , m ) + u ∞ (cid:17) . (3.14)which requires a global inversion of D (cid:101) γ w , m . The same approach is taken for LI and CLI schemes, except theprovisional velocities are obtained using local inversion only, all the inter-particle interactions are treatedexplicitly and added to the explicit term, i.e., complementary velocity (cid:101) ¯ u w − mi ; modifying Eq. (3.14) for eachvesicle, we obtain (cid:101) u w , mi = D − (cid:101) γ w , mi (cid:16) G (cid:101) γ w , mi f ( (cid:101) X w , mi , (cid:101) σ w , mi , (cid:101) λ w , mi ) + (cid:101) ¯ u w − mi (cid:17) , (3.15)where (cid:101) ¯ u w − mi is computed using (cid:101) X w , m and (cid:101) u w − mj ( j (cid:54) = i ) accounting the velocity influence from othervesicles. We only need to invert the local interaction matrices D (cid:101) γ w , mi in this scheme. Let AX + = b be the linear system that is solved at each iteration of a CLI scheme (in case of the
CLI - SDC scheme, on each of the inner step of the
SDC ). A is a block diagonal matrix,with blocks A ii corresponding to the self interactions of the i th particle. All inter-particle interactions aretreated explicitly, and thus included in the right-hand side b . We write Eq. (3.3), or Eq. (3.10), in a compactform as AX + = b + Gf + c , (3.16)0 ≤ V ( γ ; t + ) ⊥ λ ≥
0, (3.17)which is a mixed Nonlinear Complementarity Problem (
NCP ), because the
STIV function V ( γ , t ) is a nonlin-ear function of position. Note that since this is the CLI scheme, G is also a block diagonal matrix. To solve15his NCP , we use a first-order linearization of the V ( γ ; t ) to obtain an LCP and iterate until the
NCP is solvedto the desired accuracy: AX (cid:63) = b + GJ T λ . (3.18)0 ≤ V ( γ ; t + k ) + J ∆ X ⊥ λ ≥
0, (3.19)where ∆ X is the update to get the new candidate solution X (cid:63) , and J denotes the Jacobian of the volume ∇ X V ( γ , t + k ) . Algorithm 1 summarizes the steps to solve Eqs. (3.16) and (3.17) as a series of linearizationsteps Eqs. (3.18) and (3.19). We discuss the details of the LCP solver separately below.In lines 1 to 6, we solve the unconstrained system AX (cid:63) = b using the solution from previous timestep. Then, the STIV s are computed to check for collision. The loop at lines 7 −
14 is the linearized contact-resolving steps. Substituting Eq. (3.18) into Eq. (3.19), and using the fact that ∆ X = A − GJ T λ we cast theproblem in the standard LCP form 0 ≤ V + B λ ⊥ λ ≥
0, (3.20)where B = JA − GJ T . The LCP solver is called on line 9 to obtain the magnitude of the constraint force,which is in turn used to update obtain new candidate positions that may or may not satisfy the constraints.Line 13 checks the minimal-separation constraints for the candidate solution. In line 11, the collision forceis incorporated into the right-hand-side b for self interaction in the next LCP iteration. In line 14, the contactforce is updated, which will be used to form the right-hand-side b for the global interaction in the next timestep.The LCP matrix B is an M by M matrix, where M is the number of contact volumes, M = (cid:79) ( N v + N p ) .Each entry B k , p is induced change in the k th contact volume by the p th contact force. Matrix B is sparse andtypically diagonally dominant, since most STIV volumes are spatially separated.
In the contact-resolving iterations we solve an
LCP , Eq. (3.20). Most common algorithms (e.g., Lemke’salgorithm [ Lem65 ] and splitting based iterative algorithms [ Ahn81, Man77 ] ) requires explicitly formed LCP matrix B , which can be prohibitively expensive when there are many collisions. We use the minimum-mapNewton method, which we modify to require matrix-vector evaluation only, as we can perform it withoutexplicitly forming the matrix. Algorithm 1: C ONTACT - FREE TIME - STEPPING . input : X , b output : X + , f + c A ← A ( X ) b ← b ( X ) f + c ← k ← X (cid:63) ← A − b V ← getContactVolume( X (cid:63) ) while V < do J ← getContactVolumeJacobi( X (cid:63) ) λ ← lcpSolver( V ) k ← k + b ← b + GJ T λ X (cid:63) ← A − b V ← getContactVolume( X (cid:63) ) f + c ← f + c + J T λ X + ← X (cid:63)
16e briefly summarize the minimum-map Newton method. Let y = V + B λ . Using the minimum mapreformulation we can convert the LCP to a root-finding problem H ( λ ) ≡ h ( λ , y ) · · · h ( λ M , y M ) =
0, (3.21)where h ( λ i , y i ) = min ( λ i , y i ) . This problem is solved by Newton’s method (Alg. 2). In the algorithm, P (cid:65) and P (cid:70) are selection matrices: P (cid:65) λ selects the rows of λ whose indices are in set (cid:65) and zeros out all theother rows. While function H is not smooth, it is Lipschitz and directionally differentiable, and its so-calledB-derivative P (cid:65) B + P (cid:70) can be formed to find the descent direction for Newton’s method [ Erl13 ] . The matrix P (cid:65) B + P (cid:70) is a sparse matrix, and we use GMRES to solve this linear system. Since B is sparse and diagonallydominant, in practice the linear system is solved in few GMRES iterations and the Newton solver convergesquadratically.
We estimate the complexity of a single time step as a function of the number of points on each vesicle N ,number of vesicles N v . Let C N denote the cost of inverting a local linear system for one particle; then thecomplexity of inverting linear systems for all particles is (cid:79) (cid:0) C N N v (cid:1) . In [ Rah +
10, Vee + ] it is shown thatfor LI scheme C N = (cid:79) (cid:0) N log N (cid:1) . The cost of evaluating the inter-particle interactions at the N v N discretepoints using FMM is (cid:79) ( N v N ) .We assume that for each contact resolving step, the number of contact volumes is M . Assuming thatminimum map Newton method takes K steps to converge, the cost of solving the LCP is (cid:79) (cid:0) K C N N v (cid:1) ,because inverting A is the costliest step in applying the LCP matrix. The total cost of solving the
NCP problemis (cid:79) (cid:0) K K C N N v (cid:1) , where K are the number of contact resolving iterations. In the numerical simulations weobserve that the minimum map Newton method converges in a few iterations ( K ≈
15) and the number ofcontact resolving iterations is also small and independent of the problem size ( K ≈
4. Results
In this section, we present results characterizing the accuracy, robustness, and efficiency of a locally-implicittime stepping scheme combined with our contact resolution framework in comparison to schemes with nocontact resolution.
Algorithm 2: M INIMUM M AP LCP S OLVER . require : applyLCPMatrix () , V and ε output : λ e ← ε λ ← while e > ε do y ← V + applyLCPMatrix( λ ) (cid:65) ← (cid:8) i (cid:12)(cid:12) y i < λ i (cid:9) // index of active constraints (cid:70) ← (cid:8) i (cid:12)(cid:12) y i ≥ λ i (cid:9) Iteratively solve (cid:20) B − IP (cid:70) P (cid:65) (cid:21) (cid:20) ∆ λ ∆ y (cid:21) = (cid:20) − P (cid:65) y − P (cid:70) λ (cid:21) // B applied byapplyLCPMatrix τ ← projectLineSearch( ∆ λ ) λ ← λ + τ ∆ λ e ← (cid:107) H ( λ ) (cid:107) First, to demonstrate the robustness of our scheme in maintaining the prescribed minimal separa-tion distance with different viscosity contrast ν , we consider two vesicles in an extensional flow,Section 4.1. • In Section 4.2, we explore the effect of minimal separation d m and its effect on collision displacementin shear flow. We demonstrate that the collision scheme has a minimal effect on the shear displace-ment. • We compare the cost of our scheme with the unconstrained system using a simple sedimentationexample in Section 4.3. While the per-step cost of the unconstrained system is marginally lower, it re-quires a much spatial and temporal resolution in order to maintain a valid contact-free configuration,making the overall cost prohibitive. • We report the convergence behavior of different time-stepping in Section 4.4 and show that ourscheme achieves second order convergence rate with
SDC2 . • We illustrate the efficiency and robustness of our algorithm with two examples: 100 sedimentingvesicles in a container and a flow with multiple vesicles and rigid particles within a constricted tubein Section 4.5.Our experiments support the general observation that when vesicles become close, the LI scheme doesa very poor job in handling of vesicles’ interaction [ Rah + ] and the time stepping becomes unstable.The GI scheme stays stable, but the iterative solver requires more and more iterations to reach the desiredtolerance, which in turn implies higher computational cost for each time step. Therefore, the GI schemeperformance degrades to the point of not being feasible due to the intersection of vesicles and the LI schemefails due to intersection or the time-step instability. To demonstrate the robustness of our collision resolution framework, we consider two vesicles placed sym-metrically with respect to the y axis in the extensional flow u = [ − x , y ] . The vesicles have reduced areaof 0.9 and we use a first-order time stepping with LI , CLI , and GI schemes for the experiments in this test.We run the experiments with different time step size and viscosity contrast and report the minimal distancebetween vesicles as well as the final error in vesicle perimeter, which should be kept constant due localinextensibility. Snapshots of the vesicle configuration for two of the time-stepping schemes are shown inFig. 3.In Fig. 4(a), we plot the minimal distance between two vesicles over time. The vesicles continue toget closer in the GI scheme. However, the CLI scheme maintains the desired minimum separation distancebetween two vesicles. In Fig. 4(b), we show the minimum distance between the vesicles over the courseof simulation (with time horizon T =
10) versus the viscosity contrast. As expected, we observe that theminimum distance between two vesicles decreases as the viscosity contrast is increased. Consequently, for G I t = t = t = C L I Figure 3: S NAPSHOTS OF TWO VESICLES IN EXTENSIONAL FLOW USING GI AND
CLI
SCHEMES . As the distance between twovesicles decreases the configuration loses symmetry in the GI scheme as shown in the top row. Nevertheless, as shown inthe second row, the LI with minimal-separation constraint scheme maintains the desired minimum separation distanceand two vesicles also maintain a symmetric configuration. (The viscosity contrast is 500 in this simulation). a) The minimal distance of twovesicles over time − d m t M i n i m a l d i s t a n c e GICLI (b)
The minimum distance betweentwo vesicles versus the viscositycontrast − − − − Viscosity contrast M i n i m a l d i s t a n c e Figure 4: D ISTANCE BETWEEN TWO VESICLES IN EXTENSIONAL FLOW . (a) The minimum distance between two vesiclesover time for both CLI and GI schemes. The CLI scheme easily maintains the prescribed minimal separation of h. (b)The final distance (at T = ) between two vesicles as viscosity contrast is increased in GI scheme. ν ∆ t CLI LI GI e −
01 1.17 e −
01 7.66 e −
021 0.2 9.49 e −
04 9.49 e −
04 1.03 e −
031 0.1 4.49 e −
04 4.49 e −
04 4.69 e −
041 0.05 2.23 e −
04 2.23 e −
04 2.29 e −
041 0.025 1.12 e −
04 1.12 e −
04 1.14 e −
041 0.0125 5.65 e −
05 5.65 e −
05 5.68 e − e e − − e − e e −
04 1.33 e −
04 1.22 e − e e −
05 6.38 e −
05 5.96 e − e e −
05 3.05 e −
05 2.95 e − e e −
05 1.49 e −
05 1.47 e − e e −
06 7.39 e −
05 7.32 e − ν ∆ t CLI LI GI e e − T = − e − e e − − e − e e − − e − e e −
06 1.93 e −
06 1.95 e − e e −
06 9.78 e −
07 7.33 e − e e −
07 4.87 e −
07 2.01 e − e e − T = − − e e − T = − − e e − T = − − e e − − − e e − − − e e − − − Table 2: E RROR IN THE LENGTH OF VESICLES IN EXTENSIONAL FLOW . The error in the final length of two vesicles inextensional flow with respect to viscosity contrast, timestep size, and for different schemes. The experiment’s setup isdescribed in Section 4.1 and snapshots of which are shown in Fig. 3. The cases with a “ − ” indicate that either vesicleshave intersected or the GMRES failed to converge due to ill-conditioning of the system; the latter happens in the GI scheme and high viscosity contrast. Cases with superscript indicate that the flow loses its symmetry at that time. higher viscosity contrast with both GI and LI schemes, either the configuration loses its symmetry or thetwo vesicles intersect. With minimal-separation constraint, we any desired minimum separation distancebetween vesicles is maintained, and the simulation is more robust and accurate as shown in Fig. 3 andTable 2.In Table 2, we report the final error in vesicle perimeter for different schemes with respect to viscositycontrast and timestep size. With minimal-separation constraint, we achieve similar or smaller error in lengthcompared to LI or GI methods (when these methods produce a valid result). Moreover, one can use relativelylarge time step in all flow parameters — most notably when vesicles have high viscosity contrast. In contrast,the LI scheme requires very small, often impractical, time steps to prevent instability or intersection. We consider vesicles and rigid bodies in an unbounded shear flow and explore the effects of minimal separa-tion on shear diffusivity. In the first simulation, we consider two vesicles of reduced area 0.98 (to minimize19 a) t = (b) t = (c) t = (d) t = Figure 5: S HEAR FLOW EXPERIMENT SETUP . The snapshots of two vesicles in shear flow. Initially, one vesicle is placed at [ − δ ] and the second vesicle is placed at [
0, 0 ] . (a) Centroid trajectory(vesicle) − x δ t h h h h (b) Centroid trajectory(rigid) − x δ t h h h h (c) Final centroid offset d m / h ( δ ∞ ( d m ) − δ ∞ ( )) / h vesiclerigid particle Figure 6: T HE OFFSET δ t ( d m ) BETWEEN THE CENTROIDS OF TWO VESICLES IN SHEAR FLOW . The initial offset is δ = and N = discretization points are used, implying h = , where h is the distance between two discretizationpoints along vesicle boundary. (a) and (b) show δ t ( x ) for the vesicles and circular rigid particles for different minimalseparation distance d m . (c) shows the excess displacement induced by minimal separation. When collision constraintis activate (larger d m ), the particles are in effect hydrodynamically larger and the excess displacement grows linearlywith d m . the effect of vesicles’ relative orientation on the dynamics) placed in a shear flow with shear rate χ =
2. Let δ t = | y t − y t | denote the vertical offset between the centroids of vesicles at time t . Initially, two vesiclesare placed with a relative vertical offset δ as show in Fig. 5.In Fig. 6, we report δ t and its value upon termination of the simulation when x >
8, denoted by δ ∞ ( d m ) . In Figs. 6(a) and 6(b), we plot δ t with respect to the x for different minimum separation dis-tances. Based on a high-resolution simulation (with N =
128 and 8 × smaller time step), the minimaldistance between two vesicles without contact constraint is about 2.9 h for vesicles and 2.2 h for rigid par-ticles. As the minimum separation parameter d m is decreased below this threshold, the simulations withminimal-separation constraint converge to the reference simulation without minimal-separation constraint.In Fig. 6(c), we plot the excess terminal displacement due to contact constraint, ( δ ∞ ( d ) − δ ∞ ( )) / h , asa function of the minimum separation distance. When collision constraint is activate, the particles are ineffect hydrodynamically larger and the excess displacement grows linearly with d m . To compare the performance of schemes with and without contact constraints, we first consider a smallproblem with three vesicles sedimenting in a container. We compare three first-order time stepping schemes:locally implicit ( LI ), locally implicit with collision handling ( CLI ), and globally implicit ( GI ).Snapshots of these simulations are shown in Fig. 7. For the LI scheme, the error grows rapidly when thevesicles intersect and a 64 × smaller time step is required for resolving the contact and for stability. Similarly,for the GI scheme the vesicle intersect as shown in Figs. 7(e) and 7(f) and a 4 × smaller time step is neededfor the contact to be resolved. The CLI scheme, maintains the desired minimal separation between vesicles.Although the current code is not optimized for computational efficiency, it is instructive to consider therelative amount of time spent for a single time step in each scheme. Each time step in LI scheme takes about20 a) CLI t = (b) CLI t = (c) CLI t = (d) CLI t = (e) GI t = (f) GI t = Figure 7: S EDIMENTING VESICLES . Snapshots of three sedimenting vesicles in a container using three time-steppingschemes. The vesicles have reduced area of , their viscosity contrast is , and are discretized with points.The boundary is discretized with points, the simulation time horizon T = , and the time step is . (a)–(c)The dynamics of three vesicles sedimenting with CLI scheme. (d) The contact region of (c). (e) The instance when thevesicles intersect using GI scheme. (f) The contact region of (e). CLI scheme. In contrast, a single time step with the GI scheme takes, on average, 65 seconds because the solver needs up to 240 GMRES iterations to convergewhen vesicles are very close. This excessive cost renders this scheme impractical for large problems.To demonstrate the capabilities of the
CLI scheme and to gain a qualitative understanding of the scalingof the cost as the number of intersections increases, we consider the sedimentation of 100 vesicles. As thesediment progresses the number of contact regions grows to about 70 per time step. For this simulation, weuse a lower viscosity contrast of 10 and the enclosing boundary is discretized with 512 points and the totaltime is T =
30. Snapshots of the simulation with
CLI scheme are shown in Fig. 1.We observe that with first-order time stepping, both LI and CLI schemes are unstable. Therefore, we runthis simulation with second order
SDC using LI and CLI schemes. The GI scheme is prohibitively expensiveand infeasible in this case, due to ill-conditioning and large number of GMRES iterations per time-step.The LI scheme requires at least 12000 thousand time steps to maintain the non-intersection constraintand each time step takes about 700 seconds to complete. On the other hand, the CLI scheme only need1500 time steps to complete the simulation (taking the stable time step) and each time step takes about 830seconds. We repeated this experiment with 16 discretization points on each vesicle using
CLI scheme and1500 time steps are also sufficient for this case (with final length error about 3.83 e − LI scheme requires 4 × more points on each vesicle and 8 × smaller time-step size to keepvesicles in a valid configuration compared to CLI scheme.
To investigate the accuracy of the time-stepping schemes, we consider the sedimentation of three vesicleswith reduced area 0.9 in a container as shown in Fig. 7. We test LI and CLI schemes for a range of timesteps and spatial resolutions and report the error in the location of the center of mass of the vesicles at theend of simulation. The spatial resolution h is chosen proportional to time-step size and for the CLI schemethe minimal separation d m is proportional to h . As a reference, we use a fine-resolution simulation with GI scheme. The error as a function of time step size is shown in Fig. 8. As a stress test for particle-boundary and vesicle-boundary interaction, we study the flow with 25 vesiclesof reduced area 0.9, mixed with 25 circular rigid particles in a constricted tube (Fig. 9). It is well known210 − − − − Number of time steps E rr o r i n c e n t e r o f m a ss Backward Euler
CLI
Backward Euler LI first order SDC2 CLISDC2 LI second order
Figure 8: C ONVERGENCE RATE . We compare the final center of mass of three sedimenting vesicles in a container forBackward Euler and
SDC2 (second order Spectral Deferred Correction scheme). The time horizon is set to T = . Wechoose the spatial resolution proportional to the time-step. As a reference, we use the results for the GI scheme with timestep ∆ t = e − and N = discretization points on each vesicle and discretization points on the boundary. that rigid bodies can become arbitrarily close in the Stokes flow, e.g. [ LH97 ] , and without proper collisionhandling, the required temporal / spatial resolution would be very high.In this example, the vesicles and rigid particles are placed at the right hand side of the constricted tube.We use backward Euler method and search for the stable time step for schemes LI and CLI . Similarly to thesedimentation example, we do not consider the GI scheme due to its excessive computational cost.The largest stable time step for the CLI scheme is ∆ t = LI scheme, we tested the cases withup to 32 × smaller time step size, but we were unable to avoid contact and intersection between vesiclesand rigid particles.Figure 10 shows the error and minimal distance between vesicles, particles, and boundary for CLI and LI schemes with different time step sizes. Without the minimal-separation constraint, the solution divergeswhen the particles cross the domain boundary.To validate our estimates for the errors due to using a piecewise-linear approximation in the minimalseparation constraint instead of high-order shapes, we plot the minimum distance at each step for thepiecewise-linear approximation and upsampled shapes in Fig. 11. The target minimal separation distanceis set to d m = h . We observe that the actual minimal distance for smooth curve is smaller than the minimaldistance for piecewise-linear curve, while the difference between two distances is small compared to thetarget minimum separation distance.Note that the due to higher shear rate in the constricted area, the stable time step is dictated by thedynamics in that area. For such flows, we expect a combination of adaptive time stepping [ QB16 ] and thescheme outlined in this paper provide the highest speedup.
5. Conclusion
In this paper we introduced a new scheme for efficient simulation of dense suspensions of deformable andrigid particles immersed in Stokesian fluid. We demonstrated through numerical experiments that thisscheme is orders of magnitude faster than the alternatives and can achieve high order temporal accuracy.We are working on extending this approach to three dimensions and using adaptive time stepping.We extend our thanks to George Biros, David Harmon, Ehssan Nazockdast, Bryan Quaife, Michael Shel-ley, and Etienne Vouga for stimulating conversations about various aspects of this work.22 a) t = (b) t = (c) t = (d) t = Figure 9: S TENOSIS WITH VESICLES AND RIGID PARTICLES . Snapshots of vesicles and rigid particles passing aconstricted tube, the fluid flows from left to right (using proper Dirichlet boundary condition on the wall). The coloredtracers are for visualization purposes and do not induce any flow. The simulation time horizon is set to T = ,each vesicle or rigid particle is discretized with points, and the wall is discretized with points. (a) The initialconfiguration. (b)–(d) The interaction and collision between vesicles, rigid particles, and the domain boundary atdifferent instances. Without minimal-separation constraint, vesicles and rigid particles easily contact at the entranceof the constricted area. References [ Ahn81 ] B. Ahn. “Solution of nonsymmetric linear complementarity problems by iterative methods”. In:
Journal ofOptimization Theory and Applications [ All + ] J. Allard, F. Faure, H. Courtecuisse, F. Falipou, C. Duriez, and P. G. Kry. “Volume contact constraints atarbitrary resolution”. In:
ACM SIGGRAPH 2010 papers on - SIGGRAPH ’10
C (2010), p. 1. [ Alp99 ] B. K. Alpert. “Hybrid Gauss-trapezoidal quadrature rules”. In:
SIAM Journal on Scientific Computing [ Asc + ] U. M. Ascher, S. J. Ruuth, and B. T. Wetton. “Implicit-Explicit methods for time-dependent partial differ-ential equations”. In:
SIAM Journal on Numerical Analysis [ BW98 ] D. Baraff and A. Witkin. “Large steps in cloth simulation”. In:
Proceedings of the 25th annual conference onComputer graphics and interactive techniques - SIGGRAPH ’98 (1998), pp. 43–54. [ Bas + ] H. Basu, A. K. Dharmadhikari, J. A. Dharmadhikari, S. Sharma, and D. Mathur. “Tank treading of opticallytrapped red blood cells in shear flow”. In:
Biophysical Journal
23 0.5 1 1.5 2 2.5 3 3.5 410 − − − − − t M i n i m u m p a r t i c l e d i s t a n c e CLI ( ∆ t = × − ) LI ( ∆ t = × − ) LI ( ∆ t = × − ) LI ( ∆ t = × − ) − − − t E rr o r Figure 10:
ERROR FOR STENOSIS EXAMPLE . The minimal distance and error for different schemes and time step sizes. Forthe LI scheme, independently of the time step size, vesicles or particles intersect with the domain boundary, resultingin divergence. The
CLI scheme is stable and maintains the desired minimal distance. h t M i n i m a l d i s t a n c e N = Figure 11: E FFECTIVE MINIMAL DISTANCE . The minimal distance between 32-point piecewise-linear approximations(enforced by the constraint) and the true minimal distance between high-order shapes. We use an upsampled linearapproximation with N = as the surrogate for the high-order curves. [ BK06 ] A. R. Bausch and K. Kroy. “A bottom-up approach to cell mechanics”. In:
Nature Physics [ Bib + ] T. Biben, K. Kassner, and C. Misbah. “Phase-field approach to three-dimensional vesicle dynamics”. In:
Physical Review E [ Cot + ] R. W. Cottle, J.-S. Pang, and R. E. Stone.
The linear complementarity problem . Vol. 60. SIAM, 2009. [ DZ08 ] Q. Du and J. Zhang. “Adaptive finite element method for a phase field bending elasticity model of vesiclemembrane deformations”. In:
SIAM Journal on Scientific Computing [ Dur + ] C. Duriez, F. Dubois, A. Kheddar, and C. Andriot. “Realistic haptic rendering of interacting deformableobjects in virtual environments.” In:
IEEE transactions on visualization and computer graphics [ Eck + ] C. Eck, O Steinbach, and W. Wendland. “A symmetric boundary element method for contact problems withfriction”. In:
Mathematics and Computers in Simulation [ Erl13 ] K. Erleben. “Numerical methods for linear complementarity problems in physics-based animation”. In:
ACM SIGGRAPH 2013 Courses
February (2013). [ Far + ] A. Farutin, T. Biben, and C. Misbah. “3D Numerical simulations of vesicle and inextensible capsule dynam-ics”. In: (2013). [ Fau + ] F. Faure, S. Barbier, J. Allard, and F. Falipou. “Image-based collision detection and response betweenarbitrary volume objects”. In:
Proceedings of the 2008 ACM . Vol. i. 2008, pp. 155–162. [ FW05 ] K. Fischer and P Wriggers. “Frictionless 2D contact formulations for finite deformations based on themortar method”. In:
Computational Mechanics Fre07 ] J. Freund. “Leukocyte margination in a model microvessel”. In:
Physics of Fluids (1994-present) [ Fro + ] J. M. Frostad, J. Walter, and L. G. Leal. “A scaling relation for the capillary-pressure driven drainage ofthin films”. In:
Physics of Fluids [ Ghi + ] G. Ghigliotti, A. Rahimian, G. Biros, and C. Misbah. “Vesicle migration and spatial organization driven byflow line curvature”. In:
Physical Review Letters [ Gun04 ] H Gun. “Boundary element analysis of 3-D elasto-plastic contact problems with friction”. In:
Computers &structures [ Har + ] D. Harmon, E. Vouga, R. Tamstorf, and E. Grinspun. “Robust treatment of simultaneous collisions”. In:
ACM SIGGRAPH 2008 papers on - SIGGRAPH ’08 (2008), p. 1. [ Har + ] D. Harmon, E. Vouga, B. Smith, R. Tamstorf, and E. Grinspun. “Asynchronous contact mechanics”. In:
ACMTransactions on Graphics [ Har + ] D. Harmon, D. Panozzo, O. Sorkine, and D. Zorin. “Interference-aware geometric modeling”. In:
ACMTransactions on Graphics [ JJ87 ] K. L. Johnson and K. L. Johnson.
Contact mechanics . Cambridge university press, 1987. [ KL10 ] Y. Kim and M.-C. Lai. “Simulating the dynamics of inextensible vesicles by the penalty immersed boundarymethod”. In:
Journal of Computational Physics [ KW02 ] R. H. Krause and B. I. Wohlmuth. “A Dirichlet–Neumann type algorithm for contact problems with fric-tion”. In:
Computing and visualization in science [ Laa + ] A. Laadhari, P. Saramito, and C. Misbah. “Computing the dynamics of biomembranes by combining con-servative level set and adaptive finite element methods”. In:
Journal of Computational Physics
263 (2014),pp. 328–352. [ Lem65 ] C. E. Lemke. “Bimatrix equilibrium points and mathematical programming”. In:
Management science [ LH97 ] M. Loewenberg and E. J. Hinch. “Collision of two deformable drops in shear flow”. In:
Journal of FluidMechanics
338 (May 1997), pp. 299–315. [ Loe98 ] M. Loewenberg. “Numerical simulation of concentrated emulsion flows”. In:
Journal of Fluids Engineering [ Man77 ] O. Mangasarian. “Solution of symmetric linear complementarity problems by iterative methods”. In:
Jour-nal of Optimization Theory and Applications [ Mis06 ] C. Misbah. “Vacillating breathing and tumbling of vesicles under shear flow”. In:
Physical Review Letters [ Naz + ] E. Nazockdast, A. Rahimian, D. Zorin, and M. Shelley. “Fast and high-order methods for simulating fibersuspensions applied to cellular mechanics”. preprint. 2015. [ Nem + ] M. Nemer, X. Chen, D. Papadopoulos, J. Bławzdziewicz, and M. Loewenberg. “Hindered and enhancedcoalescence of drops in stokes flows”. In:
Physical Review Letters [ NW06 ] J. Nocedal and S. Wright.
Numerical optimization . Springer Science & Business Media, 2006. [ NG05 ] H. Noguchi and G. Gompper. “Vesicle dynamics in shear and capillary flows”. In:
Journal of Physics: Con-densed Matter [ OT14 ] R. Ojala and A.-K. Tornberg. “An accurate integral equation method for simulating multi-phase Stokesflow”. In: (Apr. 2014), pp. 1–22. [ Ota + ] M. a. Otaduy, R. Tamstorf, D. Steinemann, and M. Gross. “Implicit contact handling for deformable ob-jects”. In:
Computer Graphics Forum [ Pow93 ] H Power. “The completed double layer boundary integral equation method for two-dimensional Stokesflow”. In:
IMA Journal of Applied Mathematics (1993). [ PM87 ] H. Power and G. Miranda. “Second kind integral equation formulation of Stokes’ flows past a particle ofarbitrary shape”. In:
SIAM Journal on Applied Mathematics [ Poz90 ] C. Pozrikidis. “The axisymmetric deformation of a red blood cell in uniaxial straining Stokes flow”. In:
Journal of Fluid Mechanics
216 (1990), pp. 231–254. [ Poz92 ] C. Pozrikidis.
Boundary integral and singularity methods for linearized viscous flow . Cambridge Texts inApplied Mathematics. Cambridge University Press, Cambridge, 1992. [ Poz01 ] C Pozrikidis. “Dynamic simulation of the flow of suspensions of two-dimensional particles with arbitraryshapes”. In:
Engineering Analysis with Boundary Elements Pro97 ] X. Provot. “Collision and self-collision handling in cloth model dedicated to design garments”. In:
ComputerAnimation and Simulation’97 (1997). [ Pus04 ] M. A. Puso. “A 3D mortar method for solid mechanics”. In:
International Journal for Numerical Methods inEngineering [ Pus + ] M. A. Puso, T. Laursen, and J. Solberg. “A segment-to-segment mortar contact method for quadratic ele-ments and large deformations”. In:
Computer Methods in Applied Mechanics and Engineering [ QB14 ] B. Quaife and G. Biros. “High-volume fraction simulations of two-dimensional vesicle suspensions”. In:
Journal of Computational Physics
274 (2014), pp. 245–267. [ QB15 ] B. Quaife and G. Biros. “High-order adaptive time stepping for vesicle suspensions with viscosity contrast”.In:
Procedia IUTAM
16 (2015), pp. 89–98. [ QB16 ] B. Quaife and G. Biros. “Adaptive time stepping for vesicle suspensions”. In:
Journal of ComputationalPhysics
306 (2016), pp. 478–499. [ Rah + ] A. Rahimian, S. K. Veerapaneni, and G. Biros. “Dynamic simulation of locally inextensible vesicles sus-pended in an arbitrary two-dimensional domain, a boundary integral method”. In:
Journal of Computa-tional Physics [ Rah + ] A. Rahimian, S. K. Veerapaneni, D. Zorin, and G. Biros. “Boundary integral method for the flow of vesicleswith viscosity contrast in three dimensions”. In:
Journal of Computational Physics
298 (2015), pp. 766–786. [ RA78 ] J. Rallison and A. Acrivos. “A numerical study of the deformation and burst of a viscous drop in anextensional flow”. In:
Journal of Fluid Mechanics
89 (1978), pp. 191–200. [ Sac96 ] E. Sackmann. “Supported membranes: Scientific and practical applications”. In:
Science
271 (1996), pp. 43–48. [ SM94 ] A. S. Sangani and G. Mo. “Inclusion of lubrication forces in dynamic simulations”. In:
Physics of Fluids [ Soh + ] J. S. Sohn, Y.-H. Tseng, S. Li, A. Voigt, and J. S. Lowengrub. “Dynamics of multicomponent vesicles in aviscous fluid”. In:
Journal of Computational Physics [ SS01 ] S Sukumaran and U. Seifert. “Influence of shear flow on vesicles near a wall: a numerical study”. In:
Physical Review E [ Tur + ] M Tur, F. Fuenmayor, and P Wriggers. “A mortar-based frictional contact formulation for large deforma-tions using Lagrange multipliers”. In:
Computer Methods in Applied Mechanics and Engineering [ Vee + ] S. K. Veerapaneni, D. Gueyffier, D. Zorin, and G. Biros. “A boundary integral method for simulating thedynamics of inextensible vesicles suspended in a viscous fluid in 2D”. In:
Journal of Computational Physics [ Wri95 ] P. Wriggers. “Finite element algorithms for contact problems”. In:
Archives of Computational Methods inEngineering [ Wri06 ] P. Wriggers.
Computational contact mechanics . Springer Science & Business Media, 2006. [ YB12 ] A. Yazdani and P. Bagchi. “Three-dimensional numerical simulation of vesicle dynamics using a front-tracking method”. In:
Physical Review E [ YA75 ] G. K. Youngren and A. Acrivos. “Stokes flow past a particle of arbitrary shape: a numerical method ofsolution”. In:
Journal of Fluid Mechanics
69 (1975), pp. 377–403. [ ZS11 ] H. Zhao and E. S. G. Shaqfeh. “The dynamics of a vesicle in simple shear flow”. In:
Journal of FluidMechanics
674 (Mar. 2011), pp. 578–604. [ ZS13 ] H. Zhao and E. S. G. Shaqfeh. “The dynamics of a non-dilute vesicle suspension in a simple shear flow”.In:
Journal of Fluid Mechanics
725 (May 2013), pp. 709–731. [ Zha + ] H. Zhao, A. H. Isfahani, L. N. Olson, and J. B. Freund. “A spectral boundary integral method for flowingblood cells”. In:
Journal of Computational Physics [ ZP93 ] H. Zhou and C. Pozrikidis. “The flow of ordered and random suspensions of two-dimensional drops in achannel”. In:
Journal of Fluid Mechanics
255 (1993), pp. 103–127. [ ZD06 ] A. Z. Zinchenko and R. H. Davis. “A boundary-integral study of a drop squeezing through interparticleconstrictions”. In:
Journal of Fluid Mechanics
564 (Sept. 2006), p. 227. [ Zin + ] A. Z. Zinchenko, M. M. A. Rother, and R. R. H. Davis. “A novel boundary-integral algorithm for viscousinteraction of deformable drops”. In:
Physics of Fluids9.6 (1997), p. 1493.