Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method
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 , 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 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)
- et al.
Rheology and dynamics of vesicle suspension in comparison with droplet emulsion
Journal of Non-Newtonian Fluid Mechanics
(2008) An efficient numerical method for studying interfacial motion in two-dimensional creeping flows
Journal of Computational Physics
(2001)Interfacial dynamics for stokes flow
Journal of Computational Physics
(2001)- et al.
Simulating the dynamics and interactions of flexible fibers in stokes flows
Journal of Computational Physics
(2004) - et al.
A fast algorithm for hydrodynamic interactions of inextensible vesicles in 2D
Journal of Computational Physics
(2009) - et al.
Micro-macro link in rheology of erythrocyte and vesicle suspensions
Biophysical Journal
(2008) - et al.
A high-order 3D boundary integral equation solver for elliptic PDEs in smooth domains
Journal of Computational Physics
(2006) - et al.
A kernel-independent adaptive fast multipole method in two and three dimensions
Journal of Computational Physics
(2004) Hybrid gauss-trapezoidal quadrature rules
SIAM Journal on Scientific Computing
(1999)- et al.
Implicit-explicit methods for time dependent partial differential equations
SIAM Journal on Numerical Analysis
(1995)
An Introduction in Fluid Dynamics
The stress system in a suspension of force-free particles
Journal of Fluid Mechanics
Steady to unsteady dynamics of a vesicle in a flow
Physical Review Letter E
Self-diffusion of particles in shear flow of a suspension
Journal of Fluid Mechanics
Leukocyte margination in a model microvessel
Physics of Fluids
The rheological properties of suspensions of rigid particles
AIChE Journal
Transition to tumbling and two regimes of tumbling motion of a vesicle in shear flow
Physical Review Letters
Lateral migration of a two-dimensional vesicle in unbounded Poiseuille flow
Physical Review E
Integral equations of the second kind for stokes flow: direct solution for physical variables and removal of inherent accuracy limitations
Chemical Engineering Communications
Cited by (76)
TLBfind: a Thermal Lattice Boltzmann code for concentrated emulsions with FINite-size Droplets
2022, Computer Physics CommunicationsHydrodynamics simulation of red blood cells: Employing a penalty method with double jump composition of lower order time integrator
2023, Mathematical Methods in the Applied SciencesDynamics and rheology of vesicles under confined Poiseuille flow
2023, Soft MatterEffects of tunable hydrophobicity on the collective hydrodynamics of Janus particles under flows
2023, Physical Review Fluids