Adaptive Time Stepping for Vesicle Suspensions
Abstract
We present an adaptive arbitraryorder accurate timestepping numerical scheme for the flow of vesicles suspended in Stokesian fluids. Our scheme can be summarized as an approximate implicit spectral deferred correction (SDC) method. Applying a textbook fully implicit SDC scheme to vesicle flows is prohibitively expensive. For this reason we introduce several approximations. Our scheme is based on a semiimplicit linearized loworder time stepping method. (Our discretization is spectrally accurate in space.) We also use invariant properties of vesicle flows, constant area and boundary length in two dimensions, to reduce the computational cost of error estimation for adaptive time stepping. We present results in two dimensions for singlevesicle flows, constricted geometry flows, converging flows, and flows in a Couette apparatus. We experimentally demonstrate that the proposed scheme enables automatic selection of the step size and highorder accuracy.
keywords:
Vesicle suspensions, Adaptive time stepping, Spectral deferred correction, Error control, Boundary integral equations, Integrodifferentialalgebraic equations, Fluidstructure interaction, particulate flows, Stokesian flows, Moving boundary problems1 Introduction
Vesicles are deformable capsules filled with a viscous fluid. Vesicle flows refer to flow of vesicles that are suspended in a Stokesian fluid. Vesicles are differentiated from other capsules in the balance of forces in their interface. In particular, their boundary is locally inextensible (in 3D is locally incompressible). Vesicle flows find many applications many biological applications such as the simulation of biomembranes sac1996 and red blood cells ghi:rah:bir:mis2011; kao:tah:bib:ezz:ben:bir:mis2011; mis2006; nog:gom2005; poz1990.
The dynamics of a vesicle is governed by bending, tension (enforces inextensibility), hydrodynamics forces from other vesicles in the suspension, and possibly hydrodynamic forces from flow confinement boundary walls. Vesicle suspensions are modeled using the Stokes equations, a jump in the stress across each vesicle to match the interfacial forces on the membrane the vesicle, and a noslip boundary condition on the vesicles and the confinement walls.
A significant challenge in simulating vesicle flows is that their governing equations are stiff. One stiffness source is associated with the interfacial forces, bending forces, and inextensibility related tension forces (depends mainly on the curvature of the vesicle). Another stiffness source is the hydrodynamic interaction between vesicles (depends on the minimum distance between vesicles). A third stiffness source is the hydrodynamic interaction between vesicles and external boundary walls (depends on the minimum distance between the vesicles and the wall). Therefore, to resolve the dynamics of a vesicle, it may be necessary to take a small time step when the vesicle has regions of large curvature, or when a vesicle approaches a boundary wall, or approaches another vesicle. However, a large time step may be taken when a vesicle has a smooth boundary and is separated from all other vesicles and the boundary walls. These considerations necessitate that we use an adaptive timestepping scheme mostly for robustness of the code and to remove the need to manually select a time step size.
Another challenge in simulating vesicle flows is that maintaining good accuracy for long horizons requires a great number of time steps when using a loworder timestepping method. Higherorder methods in time can mitigate error accumulation over long simulation times. However, implicit higherorder methods are quite challenging to combine with adaptive schemes. The most common methods are implicit multistep schemes, but those are problematic especially for dynamics that involve moving interfaces.
In summary, the design goals for a timestepping scheme for vesicle flows is to address stiffness with reasonable computational costs, allow for adaptivity, and enable highorder time marching (at least secondorder). In the past (in other groups as well as our group), methodologies have addressed parts of the design goals above, but no method, to our knowledge, address all of them. In this paper, we propose a scheme that provides this capability.
Summary of the method and contributions
To address stiffness, we use our recent work in which wallvesicle, vesiclevesicle, and bendingtension self interactions are treated implicitly qua:bir2013b. There we used backward difference formulas (BDFs) with fixed time step sizes. As we discussed, it is possible to extend BDF with adaptive time stepping, but, since BDFs require the solution from multiple previous time step sizes but higher order methods become quite difficult to use in practice dah:lin:nev1983. We develop a new highorder method that uses spectral deferred correction to iteratively increase the order of a firstorder time integrator. By introducing SDC, we are able to iteratively construct highorder solutions that only require one previous time step.
In addition, we propose a timestepping error estimation specific to vesicle flows. The incompressibility and inextensibility conditions require that the vesicle preserves both its enclosed area and total length. Thus, errors in area enclosed by a vesicle and errors in its perimeter can be used estimate the local truncation error. In this way we avoid forming multiple expensive solutions to estimate the local truncation error.
Our contributions are summarized below.

We propose an SDC formulation to construct highorder solutions of an integrodifferentialalgebraic equation that governs vesicle suspensions.

We propose an adaptive time stepping method that uses conserved quantities to estimate the local truncation error, therefore not requiring multiple numerical solutions to form this estimate.

We conduct numerical experiments for several different flows involving multiple vesicles and confined flows that demonstrate the behavior of our scheme.
Related work
There is a rich literature on numerical methods for Stokesian particulate flows. We only discuss some representative timestepping methods for vesicle flows and we omit details on the spatial discretization. Most methods for vesicle flows are based on integral equation formulations ram:poz1998; soh:tse:li:voi:low2010; suk:sei2001; zha:sha2013b; rah:las:vee:cha:mal:moo:sam:shr:vet:vud:zor:bir2010; rah:vee:bir2010; vee:gue:zor:bir2009; vee:rah:bir:zor2011; zha:sha2011a, but other formulations based on stencil discretizations also exist laa:sar:mis2014; bib:kas:mis2005; du:zha2007; kim:lai2010.
The timestepping schemes used for vesicle flows can be grouped in three classes of methods. The first class is fully explicit methods with a penalty formulation for the inextensibility constraint. These are relatively easy to implement but inefficient due to bending stiffness. Another class of methods treats selfinteractions using linear semiimplicit time stepping. This addresses one main source of stiffness but they are fragile. No adaptive or higher than secondorder scheme have been reported for these schemes. Finally, a third class of methods treat all the vesicle interactions using linearization and semiimplicit stepping. We are aware of two papers in this direction. In zha:sha2013b a first order backward Euler scheme is used (with no adaptivity). As we mentioned in our own work qua:bir2013b; rah:vee:zor:bir2012, we used a BDF scheme for vesicle flows that treats all interactions implicitly, but again, it is not adaptive and is accurate only up to secondorder.
Finally, we briefly review the literature on SDC methods. SDC was first introduced by Dutt et al. dut:gre:rok2000. SDC has been applied to partial differential equations jia:hua2008; spe:rup:emm:bol:kra2013, differential algebraic equations bu:hua:min2012; hua:jia:min2007, and integrodifferential equations hua:lai:xia2006, but not vesicle flows. In min2003, Minion extended SDC to implicitexplicit (IMEX) methods by studying the problem where is nonstiff and is stiff. Unfortunately, our governing equations do not exhibit an additive splitting as nonstiff and stiff terms.
Algorithmically, SDC has been accelerated using preconditioners hua:jia:min2006; hua:jia:min2007 and parallel algorithms chr:ong2010; chr:mac:ong2010; chr:hay:ong2012; ong:mel:chr2012; spe:rup:emm:bol:kra2013. In our work, the SDC iterations are formed with firstorder time methods. Higherorder methods can be used to form the iterations, and this is analyzed in chr:ong:qui2009. However, such an analysis for integrodifferentialalgebraic equations is not the focus of this work.
Limitations
The main limitations of our approach are the following.

We cannot provide a proof on the convergence order of our scheme. We will see in Section 5 that each SDC correction reduces the error significantly. However, when large numbers of SDC corrections are used, the asymptotic convergence rates are unclear. Minion observes similar behavior for stiff problems in min2003. Understanding the asymptotic rates of convergence could be used to further optimize the procedure for choosing optimal time step sizes.

We have not explored the scheme for vesicle suspensions where the viscosity inside and outside each vesicles differs (viscosity contrast). The difficulty arises because flows with a viscosity contrast include a doublelayer potential of the fluid velocity.

Currently, each vesicle must use the same time step size. Therefore, in flows with multiple vesicles, the time step size is controlled by the vesicle requiring the smallest time step. This can result in the actual global error being much less than the desired error (see Couette example). Simulations could be accelerated if each vesicle could have its own time step size.
Let us remark that we do not use adaptivity in space.
Outline of the paper
In Section 2, we briefly discuss SDC for initial value problems (IVPs), and then discuss how to extend SDC to vesicle suspensions. In Section 3, we discretize the differential and integral equations arising in the SDC framework, discuss preconditioning, and provide complexity estimates. Section 4 discusses our adaptive time stepping strategy, and numerical results are presented in Section 5. Finally, we make concluding remarks in Section 6, and we reserve a more in depth formulation of our governing equations in the SDC framework Appendix A.
Notation
In Table 1, we summarize the main notation used in this paper.
Symbol  Definition 

Boundary of vesicle  
Parameterization of at time and parameterized in arclength  
Tension of vesicle at time and parameterized in arclength  
Boundary of the confined geometry  
Bending operator due to vesicle  
Tension operator due to vesicle  
Surface divergence operator due to vesicle  
Singlelayer potential due to vesicle and evaluated on vesicle  
Doublelayer potential due to and evaluated on vesicle  
Background velocity due to a farfield condition  
Velocity of vesicle due to hydrodynamic forces from vesicle  
Residual of equation (7) due to the provisional solution  
Error between the exact vesicle position and the provisional position  
Error between the exact vesicle tension and the provisional tension  
Area and length of the vesicles at time  
Error in area and length  
Desired final tolerance for the area and length of the vesicles  
Maximum amount the time step is allowed to be scaled up per time step  
Minimum amount the time step is allowed to be scaled down per time step  
Multiplicative safety factor for the new time step size 
2 Formulation
In this section, we briefly summarize the SDC formulation for IVPs and then extend SDC to an integrodifferential equation that governs vesicle dynamics. We only discuss the theory of SDC and outline its numerical implementation in Section 3.
2.1 Spectral Deferred Correction
In its original development dut:gre:rok2000, SDC iteratively constructed a highorder solution of the IVP
(1) 
While classical deferred correction methods discretize the time derivative in (1), SDC uses a Picard integral to avoid unstable numerical differentiation. Equation (1) is reformulated as
(2) 
Any time stepping scheme can be used to generate a provisional solution . Given this provisional solution of (2), the residual is defined as
(3) 
The integral in (3) is approximated by a order accurate quadrature rule. The error satisfies
(4) 
and an approximation of is computed, generally using the same numerical method used to compute . Finally, the provisional solution is updated to and the procedure is repeated as many times as desired. The process of forming the residual , approximating the error , and updating the provisional solution is referred to as an SDC correction or iteration.
The accuracy of SDC depends on the discretization of (2) and (4) and the accuracy of the quadrature rule required to evaluate . An abstract error analysis for applying deferred correction methods to the operator equation has been preformed in boh:ste1984; lin1980; ske1982; ste1973. In han:str2011, this abstract framework was applied to (1) by defining , and the main result is state in Theorem 4.2. Here we only summarize the result to avoid introducing additional notation.
Theorem 2.1.
Let have bounded derivatives, and be the unique solution of (1). Suppose that the time integrator is stable. That is, there exists some which only depends on and , such that for all ,
Moreover, suppose that is of order meaning that , where . Suppose that a numerical solution of satisfies
where depends only on the derivatives of and on , but not on . Assuming the residual is computed exactly, if is formed with the order time integrator to approximate , then
However, since is approximated with a order quadrature, the asymptotic error is
Theorem 2.1 tells us that by estimating the error with a firstorder method, which is the only order we consider in the SDC framework, the order of accuracy is increased by one, with the constraint that this convergence is limited by the accuracy of the quadrature rule for approximating (3). However, the theorem states nothing about the stability of SDC. In dut:gre:rok2000, the authors consider three discretizations of (2) and (4): fully explicit, fully implicit, and linear combinations of the two. We do not consider explicit methods since the governing equations are stiff. The other two methods could be used for vesicle suspensions, but both would require solving nonlinear equations.
We have successfully used a variant of IMEX methods for vesicle suspensions qua:bir2013b; rah:vee:bir2010; vee:gue:zor:bir2009, and we couple these methods with SDC in this work. IMEX methods asc:ruu:wet1995 are a family of time integrators that treat some terms (generally, linear) implicitly and other terms explicitly. IMEX methods for additive splittings of (1) were first applied to SDC by Minion in min2003. We summarize the numerical results from min2003 since their behaviour resembles results that we will observe. First, the Van der Pol oscillator is considered in a nonstiff, mildly stiff, and stiff regime, and the number of SDC iterations ranges from two to six. In the nonstiff regime, the error behaves according to Theorem 2.1. In the mildly stiff case, the correct asymptotic result is observed, but not until is much smaller. In the stiff regime, the convergence behaviour differs considerably from the formal error. This behaviour is attributed to order reduction which is further analyzed. The author proceeds to claim that “the correct asymptotic convergence rates would be observed given sufficiently small ; however, this is not the relevant issue in most applications”. Two other examples considered are a system of differential equations resulting from a spatial discretization of an advectiondiffusion equation, and Burgers equation. Convergence rates up to fifthorder are achieved. As before, as the system becomes stiffer, smaller time steps are required before the expected convergence behaviour is observed.
2.2 SDC for Vesicle Suspensions
Let be a collection of vesicles parameterized by and with tension (Figure 1). We use the integrodifferential equation from vee:gue:zor:bir2009 to balance the bending and tension forces of the vesicle with the stress jump of the fluid across the vesicle membrane, and an algebraic constraint to enforce the inextensibility condition. We start by defining the velocity of vesicle due to the hydrodynamic forces from vesicle
where is the Kronecker delta function,
and is the background velocity (unconfined flows) or the velocity due to solid walls (confined flows). In the case of confined flows, we use the doublelayer potential of an unknown density function defined on the boundary of the solid walls. The extra equation comes from a nonslip boundary condition on the solid walls and the details are presented in rah:vee:bir2010. We point out that is not symmetric, meaning that for all , and , , are all linear in their second argument.
The notation we are using is chosen so that terms such as , which approximates , are understood to mean
where is the derivative of with respect to its parameterization variable. The tension acts as a Lagrange multiplier to satisfy the inextensibility constraint
(5) 
where
which is also linear in its second argument. Equation (5) can be eliminated using the Schur complement of the tension to write in terms of the positions , . Then, can be written entirely in terms of and , and the noslip boundary condition of vesicle gives
(6) 
The formulation (6) is easiest to present our numerical methods, but we in fact do not eliminate the tension in our implementation. The resulting changes to the SDC formulation are presented in Appendix A.


Following the SDC method we reformulate (6) as
(7) 
We form a provisional solution at a set of quadrature nodes in using a method that we describe in the next section (Section LABEL:s:vesiclePicard), and then evaluate
(8) 
with a quadrature rule that we also define in the next section (Section LABEL:s:vesicleQuadrature). Then, the error satisfies
which we write using the residual as
(9) 
We define the new provisional solution as , where is an approximate solution of (9). Again, this procedure of computing the residual (8), numerically solving (9) for the error, and updating the provisional solution is what we call an SDC iteration. Assuming that we take firstorder SDC iterations and is the quadrature error for computing the residual , from Theorem 2.1, we expect that the asymptotic rate of convergence is .
3 Numerical Scheme
4 Adaptive Time Stepping
5 Results
6 Conclusions
Appendix A Complete SDC Formulation
\biboptionssort&compress