Elsevier

Journal of Computational Physics

Volume 229, Issue 18, 1 September 2010, Pages 6466-6484
Journal of Computational Physics

Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method

https://doi.org/10.1016/j.jcp.2010.05.006Get rights and content

Abstract

We consider numerical algorithms for the simulation of hydrodynamics of two-dimensional vesicles suspended in a viscous Stokesian fluid. The motion of vesicles is governed by the interplay between hydrodynamic and elastic forces. Continuum models of vesicles use a two-phase fluid system with interfacial forces that include tension (to maintain local “surface” inextensibility) and bending. Vesicle flows are challenging to simulate. On the one hand, explicit time-stepping schemes suffer from a severe stability constraint due to the stiffness related to high-order spatial derivatives in the bending term. On the other hand, implicit time-stepping schemes can be expensive because they require the solution of a set of nonlinear equations at each time step.

Our method is an extension of the work of Veerapaneni et al. [S.K. Veerapaneni, D. Gueyffier, D. Zorin, G. Biros, A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D, Journal of Computational Physics 228(7) (2009) 2334–2353], in which a semi-implicit time-marching scheme based on a boundary integral formulation of the Stokes problem for vesicles in an unbounded medium was proposed.

In this paper, we consider two important generalizations: (i) confined flows within arbitrary-shaped stationary/moving geometries; and (ii) flows in which the interior (to the vesicle) and exterior fluids have different viscosity. In the rest of the paper, we will refer to this as the “viscosity contrast”. These two problems require solving additional integral equations and cause nontrivial modifications to the previous numerical scheme. Our method does not have severe time-step stability constraints and its computational cost-per-time-step is comparable to that of an explicit scheme. The discretization is pseudo-spectral in space, and multistep BDF in time. We conduct numerical experiments to investigate the stability, accuracy and the computational cost of the algorithm. Overall, our method achieves several orders of magnitude speed-up compared to standard explicit schemes.

As a preliminary validation of our scheme, we study the dependence of the inclination angle of a single vesicle in shear flow on the viscosity contrast and the reduced area of the vesicle, the lateral migration of vesicles in shear flow, the dispersion of two vesicles, and the effective viscosity of a dilute suspension of vesicles.

Introduction

Vesicles are closed lipid membranes suspended in a viscous medium. The mechanical deformation of vesicles and their interaction with viscous fluids are thought to play an important role in many biological phenomena [13], [31] and are used experimentally to understand properties of biomembranes [30]. In addition, vesicle mechanics have been used as models for the motion of red and white blood cells [19], [22], whose quantitative description will help in better understanding blood rheology.

In this article, we focus on numerical schemes for continuum models of vesicle dynamics. This is a challenging problem because the motion and shape of the vesicles must be determined dynamically from a balance of interfacial forces with fluid stresses. The shape dynamics of fluid vesicles is governed by the coupling of the flow within the membrane of the vesicle with the hydrodynamics of the surrounding bulk fluid. Following our previous work on vesicle flows [34], we present a semi-implicit numerical scheme for the simulation of the motion of arbitrarily shaped vesicles that can have a viscosity contrast with the background fluid. We also extend our formulation to handle interior flows and interaction of vesicles with other moving particles with prescribed motion.

Our method is based on an integral equation formulation. In particulate flow problems involving vesicles, the elastic and incompressibility properties of their membranes must be resolved and the numerical schemes must be modified in order to accommodate these properties and to solve the resulting set of equations. Details of the boundary integral formulation for elastic interfaces and incompressible vesicles can be found in the works of Pozrikidis [24], [25].

The overwhelming majority of works on particulate flows uses explicit schemes that pose severe restrictions on the time step. In contrast, semi-implicit methods result in two to three orders of magnitude larger step size that is almost independent of the spacial grid size [34]. In contrast to stencil-based methods (e.g. finite element methods), integral equation formulations avoid discretization of the overall domain and instead discretize only the vesicle-boundary and the boundary of the enclosing domain. This is the main reason that integral equations have been used extensively for vesicle, and more generally, particulate and interfacial flow simulations [25].

The boundary integral formulation coupled to the shape dynamics results in an integro-differential equation that is constrained by the local inextensibility. Extending our previous work [34], we use semi-implicit time-stepping, fast summation schemes, and spectral discretization in space. The combination of these approaches for flows with interface singularities is not unique. However, we are unaware of any previous analysis and application of implicit time-stepping schemes combined with fast solvers to vesicles that have a viscosity contrast with the surrounding fluid and are interacting with confined boundaries. These improvements enable the simulation a large number of interacting vesicles, as described in Sections 3 Numerical algorithm, 4 Numerical experiments, and depicted in Fig. 1.

The main contributions of this paper are:

  • The extension of the techniques developed in [14], [20], [34] to vesicle flows in confined geometry and vesicles with viscosity contrast.

  • The numerical investigation of the stability and accuracy of the time-stepping schemes.

  • A preliminary validation of our methodology by comparing our numerical results to results in the literature.

In particular, for validation, we investigate (i) the dependence of vesicles’ inclination angle in shear flow on viscosity contrast and reduced area; (ii) the lateral migration of vesicles in shear flow due to collision; and (iii) the rheology of a dilute suspension of vesicles.

The most significant limitation of our method is that the number of Fourier modes used to represent the vesicle membrane and the time step are not chosen adaptively. The former is a minor limitation (in 2D) but the latter is quite significant. Our spectral discretization (which we combine with a special high-order scheme for singular integrals) in space [34] results in discretization errors that are dominated by the time-stepping scheme. In our experiments, 64–128 Fourier modes in space are typically sufficient to fully resolve the shapes of the vesicles in the flow regimes we have examined. For more concentrated suspensions, adaptive schemes combined with a posteriori estimates may be necessary.

We solve the discretized system of equations using the Generalized Minimum Residual Method (GMRES) [29] with an appropriate set of preconditioners, which are based on the spectral properties of the operators. Nonetheless, for very small viscosity contrasts ν  1 (see Table 1 for its definition), the spectral properties of the operators change and a generic preconditioner, as we use here, fails to fully compensate for the poor conditioning of the operators.

Vesicles are used, theoretically and experimentally, to investigate the properties of biological membranes [30], blood cells [19], [22], and drug-carrying capsules [32].

Integral equation methods have been used extensively for the simulation of Stokesian particulate flows, mostly for droplets (with or without viscosity contrast) and bubbles. These methods were introduced by Youngren and Acrivos [38] for a flow past a rigid particle of arbitrary shape. An excellent review of numerical methods for Stokesian flows is done by Pozrikidis [25]. The present work is based on a formulation derived by Rallison and Acrivos [27] for two fluids separated by an interface with surface forces and the work of Power et al., who, in a series of papers [20], [21], introduced an integral equation formulation of the Stokes problem on a multiply-connected domain with Dirichlet boundary conditions.

In spite of the large body of literature devoted to investigating the rheological properties of red blood cell and vesicles suspensions, to the best our knowledge, the work on numerical methods for vesicle flows with viscosity contrast and confined boundaries is rather limited. Freund [8] considers vesicles with no viscosity contrast in a bounded domain. In his work, the boundaries are treated as panels fixed to their location with virtual springs. Zhou and Pozrikidis [39] consider the flow of a periodic suspension of 2D viscous drops between two parallel plane walls, for which an explicit expression of Green’s function is available.

Let us also mention works related to the test-case flows we have used to validate our numerical method. Kraus et al. [13] studied the dynamics of a vesicle and its steady-state inclination angle in the absence of viscosity contrast. Beaucourt et al. [5] tackled the same problem in the presence of viscosity contrast. Kantsler and Steinberg [10] reported results from experimental study of the inclination angle of vesicles and their transition from tank-treading to tumbling. Misbah [18] looked at the theoretical aspects of a vesicle’s inclination angle problem. Loewenberg [15], Loewenberg and Hinch [16] studied the dispersion of drops in shear flow and Eckstein et al. [7] investigated particle–particle interaction and their lateral migration in an experimental setting. Rheology of (dilute) suspension of vesicles was investigated by Danker et al. [6], Loewenberg [15], Ramanujan and Pozrikidis [28], and Vitkova et al. [35]. Danker et al. [6] investigated the rheological properties of a dilute suspension of vesicles in shear flow analytically.

We propose a computational scheme for the evolution of vesicles in a confined domain. The scheme builds upon our previous work [34]. We also extend our method to the case where there is a viscosity contrast between the suspending fluid and the internal fluid of the vesicle. Our scheme is based on Lagrangian tracking of marker particles on the vesicle, semi-implicit time discretization and spectral representation of the interface, together with high-order accurate quadrature rules. These choices result in a spectrally accurate method in space and second-order accurate method in time.

High-order accuracy in space is ensured by using a Fourier basis discretization for all functions and computing derivatives in Fourier domain, as well as high-order, Gauss-trapezoidal quadrature rules [1] for discretization of single-layer potentials. In time, we use a semi-implicit marching scheme [2]. This discretization yields a linear system of equations for each time step, which is solved using GMRES. One significant challenge in simulating vesicle dynamics is the numerical stiffness of the governing integro-differential equations [34]. To gain insight on the spectral properties of operators we use a “frozen coefficient” analysis on the unit circle. This analysis allows us to construct a preconditioner for the GMRES solver. Putting everything together, we were able to achieve high accuracy in space and time, while taking large time steps without incurring high computational costs. Our formulation for confined boundaries is based on the method in [20]. Finally, we resolve nearly-singular integrals, which arise when vesicles come close the fixed boundary, using the method proposed in [36].

Throughout this paper, lower case letters will refer to scalars, and lowercase bold letters will refer to vectors. We use ⊗ for the tensor product of two vectors and ∣·∣ to denote the measure of its argument (e.g. the Euclidean norm of a vector or the area of a domain). We denote the jump across interfaces by 〚u  u+  u, where u±(x)  limh↓0u(x ± hn), n denoting the outward normal to the boundary. We denote the convolution of an integral kernel K with density η by K[y,η](x)ΓK(x,y)η(y)ds(y), where the product inside the integral should be interpreted as a tensor operation when K is a tensor and as a dot product when K is a vector. In Table 1, we list symbols and operators frequently used in this paper.

Finally, the rest of the paper is organized as follows: In Section 2, we state the problem and its formulation. In Section 3, we outline the numerical scheme we use to solve the derived equations. In Section 4, we report results from numerical experiments we performed to demonstrate the stability of the proposed time-marching scheme in different flow regimes and geometry configurations. In particular, we investigate the effect of fixed boundaries on the time-stepping numerical stability of our scheme. We conclude in Section 4.4 with a discussion of the rheology of dilute suspensions of vesicle flows with viscosity contrast.

Section snippets

Formulation

Consider a suspension of vesicles in an ambient Newtonian fluid. Let Ω be the domain of interest, an open bounded subset of R2 that can be multiply-connected and whose boundary consists of M + 1 infinitely differentiable curves Γ0,  , ΓM, among which, Γ0 denotes the enclosing boundary of the domain. The ambient fluid has viscosity μ0. The vesicles are evolving under the influence of an imposed velocity field. Let γp(p = 1,  , N) denote the boundary of the pth vesicle, ωp denote the domain enclosed by γp

Numerical algorithm

We use a multistep time-marching scheme. In Section 3.1, we give the details of our time-stepping schemes and in Section 3.2, we outline our approach to spatial discretization. Following [34], we adopt a Lagrangian formulation, which simplifies the implementation of the high-order multistep schemes.

Numerical experiments

In this section, we present results on the convergence, stability, and algorithmic complexity of the proposed methods, which we have implemented in MATLAB. We preform the following tests:

  • We consider a single vesicle in Section 4.1. Our goal is to demonstrate the stability and accuracy of our scheme as a function of the parameters, specifically viscosity contrast ν. We also report the dependence of tank-treading inclination angle on a vesicle’s viscosity contrast and shear rate.

  • We consider the

Conclusions

We proposed a semi-implicit numerical scheme to simulate the motion of inextensible vesicles suspended in bounded or unbounded domains. For several test-cases, we have demonstrated, through the use of numerical experiments, that the proposed scheme does not exhibit a mesh-dependent high-order stability constraint on the time-step size. Our scheme exhibits second-order accuracy in time and spectral-accuracy in space. We have presented efficient low-cost preconditioners to solve the discrete

References (39)

  • G.K. Batchelor

    An Introduction in Fluid Dynamics

    (1970)
  • G.K. Batchelor

    The stress system in a suspension of force-free particles

    Journal of Fluid Mechanics

    (1970)
  • J. Beaucourt et al.

    Steady to unsteady dynamics of a vesicle in a flow

    Physical Review Letter E

    (2004)
  • E.C. Eckstein et al.

    Self-diffusion of particles in shear flow of a suspension

    Journal of Fluid Mechanics

    (1977)
  • Jonathan B. Freund

    Leukocyte margination in a model microvessel

    Physics of Fluids

    (2007)
  • D.J. Jeffrey et al.

    The rheological properties of suspensions of rigid particles

    AIChE Journal

    (1976)
  • V. Kantsler et al.

    Transition to tumbling and two regimes of tumbling motion of a vesicle in shear flow

    Physical Review Letters

    (2006)
  • B. Kaoui et al.

    Lateral migration of a two-dimensional vesicle in unbounded Poiseuille flow

    Physical Review E

    (2008)
  • S.J. Karrila et al.

    Integral equations of the second kind for stokes flow: direct solution for physical variables and removal of inherent accuracy limitations

    Chemical Engineering Communications

    (1989)
  • Cited by (76)

    View all citing articles on Scopus
    View full text