Two-dimensional lattice Boltzmann simulations of vesicles with viscosity contrast

We present a numerical approach to simulate the dynamics of viscous vesicles (their internal and external fluids have different viscosities). The flow is computed using the lattice Boltzmann method and the fluid-vesicle two-way coupling is achieved using the immersed boundary method. The viscosity contrast (defined as the ratio of the internal to the external viscosities) is included using a geometrical algorithm that detects if a fluid node is either located inside or outside a vesicle. Our two-dimensional simulations successfully reproduce the tank-treading and tumbling dynamical states known for a viscous vesicle when it is subjected to simple shear flow. A good qualitative agreement between our simulation results and literature data is obtained. Moreover, we quantitatively analyze how inertia influences the dynamics of a vesicle and as an outlook we present an application of our method to the flow of multiple viscous vesicles in a microfluidic constriction.


Introduction
A vesicle is a fluid-filled particle with a membrane made of lipid molecules. It behaves as a compartment that encapsulates a fluid (or a suspension) and protects it from an external suspending fluid. Vesicles are often used for microencapsulation of active materials for drug delivery or as a biomimetic model to study biological cells, for example, blood cells. One of the open questions that attracted the interest of scientists in recent decades is: How does a vesicle deform and behave dynamically when it is transported by flow? There is a rich literature dealing with this problem theoretically (Keller and Skalak 1982;Seifert 1999;Misbah 2006;Vlahovska and Gracia 2007;Lebedev et al. 2007;Kaoui et al. 2009;Abreu and Seifert 2013), numerically (Kraus et al. 1996;Biben and Misbah 2003;Beaucourt et al. 2004;Noguchi and Gompper 2005;Du et al. 2006;Finken et al. 2008;Veerapaneni et al. 2009;Kaoui et al. 2011;Zhao and Shaqfeh 2011;Kaoui et al. 2012;Maitre et al. 2012;Laadhari et al. 2012;Salac and Miksis 2012;Doyeux et al. 2013;Krüger et al. 2013;Halliday et al. 2013), and experimentally (De Haas et al. 1997;Mader et al. 2006;Kantsler and Steinberg 2005;2006). The present paper is an additional contribution to this field. It reviews a numerical method capable of simulating multiple vesicles with viscosity contrast (the viscosities of their internal and external fluids are different η int = η ext ) and flowing in complex geometries.
Mathematically, the problem of vesicle dynamics consists of tracking the motion of a free-moving boundary (the membrane) separating two distinct spacial domains (the internal and the external fluids). This is a challenging problem: The location of the membrane is not known a priori and its motion results from its interaction with the surrounding fluid flow. This latter itself is also unknown since it is affected by the time-dependent motion and deformation of the membrane. The problem of vesicle dynamics falls in the class of tracking interfaces and fluid-structure interaction (Krüger et al. 2013). Several numerical approaches were proposed to fully resolve this problem while considering the two-way coupling between the fluid flow and the membrane dynamics. There are methods based on front tracking (called also Lagrangian methods), such as, the immersed boundary method (IBM) combined with a fluid solver (Kaoui et al. 2011;Kaoui et al. 2012;Noguchi and Gompper 2005;Finken et al. 2008) and the boundary integral method (BIM) (Kraus et al. 1996;Beaucourt et al. 2004;Veerapaneni et al. 2009;Zhao and Shaqfeh 2011). Alternative approaches are based on front capturing (named also Eulerian methods), for example, the phase-field method (PFM) (Biben and Misbah 2003;Beaucourt et al. 2004;Du et al. 2006;Ghigliotti et al. 2010) and the level-set method (LS) (Maitre et al. 2012;Laadhari et al. 2012;Salac and Miksis 2012;Doyeux et al. 2013). Some of these methods are more adequate for simulating the situation of a vesicle immersed in an unbounded fluid, such as the BIM. Other methods, like the IBM combined with a fluid solver, are more adequate to simulate situations where a vesicle is confined and flowing in complex geometries (e.g., in the blood microcirculatory system or in microfluidic devices).
In this paper, we couple the vesicle dynamics and the fluid flow using the IBM. The flow is computed with the lattice Boltzmann method (LBM) (Succi 2001;Sukop and Thorne 2006). Other groups have already proposed numerical methods for vesicle dynamics in which the flow is computed with the LBM (e.g., Dupin et al. (2007) and Halliday et al. (2013)). Here, we present a method that allows to simulate the dynamics of vesicles with a viscosity contrast. We implement the viscosity contrast using a simple geometrical algorithm that detects if a fluid node is located either inside or outside a vesicle. Our approach is simpler, in the sense that it uses only a single phase/component fluid flow and does not require solving any additional field for the viscosity contrast in comparison to front capturing methods (Biben and Misbah 2003;Beaucourt et al. 2004;Du et al. 2006;Maitre et al. 2012;Laadhari et al. 2012;Salac and Miksis 2012;Doyeux et al. 2013;Halliday et al. 2013). Despite the simplicity of our algorithm, it has demonstrated its ability to capture the essential features known for the dynamics of a vesicle under shear flow. In our previous paper (Kaoui et al. 2012), we have already applied this approach to study the effect of wall confinement on the transition between the dynamical states of a single vesicle subjected to simple shear flow. In the current paper, we analyse in detail how the simulation data are sensitive to the parameter space. In particular, the inertia present in our approach via the LBM has an effect despite the low values of the Reynolds number we use (Re < 1). We apply our method to the case of multiple viscous vesicles flowing in a microfluidic constriction. This is a suitable application to demonstrate the strength of our method to handle at once the deformation of the particles (with the IBM), the viscosity contrast and the flow through complex geometries (thanks to the LBM).
This paper is organized as follows: First, we introduce our computational method. Second, we show that our approach is capable to capture the known dynamics of a single viscous vesicle subjected to shear flow. We briefly discuss the effect of inertia. Third, we present simulations of multiple viscous vesicles flowing through a microfluidic constriction. Finally, we conclude.

Lattice Boltzmann method (flow)
The flow of the internal and the external fluids follows the well-known Navier-Stokes (NS) equations, where ρ and η are the mass density and the dynamic viscosity of the fluid (the hydrodynamical properties of the fluid). u and p are its velocity-and pressure fields (the unknown quantities of our problem). The above equations are solved while considering boundary conditions on the vesicle membrane ∂Ω. However, the location of this is on its own an additional unknown quantity of the problem. Its motion results from its hydrodynamical interaction with its surrounding fluid flow. Therefore, to obtain the dynamics of the vesicle, we need to acquire information about the flow. In the present paper, instead of solving the above NS equations, we use an alternative mesoscopic approach: the lattice Boltzmann method (LBM), that recovers the solution of the NS equations in the limit of low Mach and Knudsen numbers. The LBM has gained popularity among scientists and engineers because of its relatively straightforward implementation, in comparison to other classical computational fluid dynamics solvers, such as, the finite element method or the finite volume method. For a review, see, e.g., Refs. (Succi 2001;Sukop and Thorne 2006). The fluid is considered as a collection of pseudo-particles that live and move on a lattice under the action of external forces. The single particle distribution function f i (r, t) gives the probability to find an elementary portion of the fluid at a position r with a velocity in the direction e i . In the LBM, both the position and the velocity space are discretized. Here, we use nine discrete velocity directions in 2D resulting in the socalled the D2Q9 lattice. Δx is the lattice spacing and e i are the discretized velocity directions. The time evolution of f i is governed by the lattice Boltzmann equation, where Δt is the time step (in the following we choose Δt = Δx = 1. Δt = 1 is the smallest time scale involved in the physical problem we study below). The lattice Boltzmann equation consists of two parts: (i) the advection part (the left-hand side of Eq. 2) and (ii) the collision part (the right-hand side of Eq. 2). Ξ i is the collision operator specifying the collision rate between the fluid pseudo-particles. It can be approximated by the Bhatnagar-Gross-Krook (BGK) operator (Succi 2001), Ξ i = −(f i − f eq i )/τ , which describes the relaxation of f i towards its local equilibrium, f eq i , on a time scale τ . This relaxation time is related to the macroscopic dynamical viscosity η via η = ρc 2 s τ − 1 2 where c s = 1/ √ 3 is the lattice speed of sound. The equilibrium distribution is given by a truncated Maxwell-Boltzmann distribution (Succi 2001), where the ω i are the weight factors resulting from the velocity space discretization. The hydrodynamical macroscopic quantities are computed using the first and the second moments of f i : (i) the local pressure p = ρc 2 s = c 2 s i f i and (ii) the local velocity u = i f i e i /ρ. F, on the right hand side of Eq. 2, is an external applied body force. In the context of the current contribution, it is the membrane force as it is detailed in the following subsection.

Vesicles (structure)
Vesicles are closed lipid membranes. They are commonly used as a biomimetic model for biological cells, for example, to study the dynamics and rheology of red blood cells. They are also used for microencapsulation of active materials for drug delivery by the pharmaceutics industry. The thickness of their membrane (∼ 5 nm) is negligibly small as compared to their size (typically ∼ 10 μm). Therefore, vesicles are modeled as a closed surface with zero thickness. Vesicles have two main properties: their volume and area are fixed, i.e., they do not change in time when the vesicle undergoes mechanical deformation. The constraint of volume conservation is usually fullfiled if the encapsulated fluid is an incompressible Newtonian fluid. The second constraint of the surface area conservation is due to the resistance of the membrane to undergo compression/extension deformation. This is achieved mathematically by using a local Lagrangian multiplier σ that plays the role of an effective surface tension. It is a dynamical local quantity that evolves in time, as such as to keep the vesicle area constant. The associated tension energy is given by E S = ∂Ω σ (s) ds, where the integration is performed over the membrane surface ∂Ω.
The vesicle membrane experiences also resistance to bending, E B = κ B 2 ∂Ω c 2 (s) ds, where κ B is the bending rigidity and c(s) the membrane local curvature. When the membrane is bent due to applied hydrodynamic stresses, it tries to regain its equilibrium configuration by exerting a reaction force back on its surrounding fluid. The total membrane force is derived by taking the functional derivative F = δ(E B + E S )/δr. In 2D, it is given by Kaoui et al. (2011) where n and t are, respectively, the unit normal and tangent vectors and s is the curvilinear coordinate. The geometrical quantities on the membrane (the local curvature c and its derivatives) required to evaluate the force (4) are computed using the finite difference method. The vesicle membrane is a freely-moving boundary whose shape and position are not known a priori. The membrane dynamics results from its interaction with the surrounding fluid flow. In the present work, the motion of the membrane is computed using a front tracking method in 2D. The membrane is represented by a moving Lagrangian mesh (a contour in 2D), and the fluid flow is computed on a fixed Eulerian lattice. Instead of solving the Navier-Stokes equations, inside and outside the vesicle, we use alternatively the LBM to compute the flow (see the previous subsection). The vesicle model can be easily adapted to other kinds of deformable particles such as droplets or capsules by using the appropriate constitutive laws instead of Eq. 4. The membrane force F acting on the fluid is included in the model by adding the term (F · e i ) to the right-hand side of Eq. 2.

Fluid-structure coupling
The vesicle membrane is modeled as a continuum medium characterized by macroscopic mechanical properties, such as the membrane rigidity κ B . The way the membrane reacts back to the external forces (e.g, hydrodynamic stresses) results from the response and the rearrangement of its microscopic constituents (the lipid molecules) that we do not cope with directly in our approach. The membrane is considered to be sharp (zero-thickness), and it is represented with a cluster of marker points interconnected via elastic springs. This constitutes a moving Lagrangian mesh (structure) for which we need to predict the motion during time.
The tension built up between the two adjacent points, separated by the distance Δs(s, t) at the time t, is given by where κ S is the spring constant and Δs(s, t = 0) is the initial referential distance between the two adjacent points. The membrane mesh is immersed in a fixed Eulerian mesh representing the fluid and its dynamics results from its action-reaction with the flow. In the present work, we follow Peskin's immersed boundary method (IBM) to accomplish the two-way coupling between the fluid flow and the membrane dynamics (Peskin 1977). The physical quantities computed in each mesh are matched together by interpolation. This method consists of two main steps: (i) Advection: First, the flow is computed on the whole fixed Eulerian mesh using the LBM. Then, the velocity v at every position r m belonging to the membrane is computed by interpolating, at this point, the velocities u of its nearest surrounding fluid nodes where D is a function suggested by Peskin (2002). It has a peak on the membrane point r m and decays at a distance equal to twice the lattice spacing (2Δx) after which it vanishes: for |x| 2Δx and |y| 2Δx. D(r) = 0 elsewhere. We choose this form of D because it is expected to increase the accuracy of the method as discussed in Krüger et al. (2011). The above interpolation can be used only under the assumption of a continuous velocity field across the membrane, (8) and considering the mesh points to be massless (as tracer particles with no effect on the flow at this stage). Finally, each membrane point is advected using an Euler scheme: (ii) Reaction: Once the membrane is advected, it gets deformed. The resulting shape of the vesicle is not necessary its equilibrium shape which would minimize its energy. Thus, the membrane tries to minimize the energy by exerting a reaction force (4) back on its surrounding fluid. This force is distributed to the fluid nodes using the following expression, Every fluid point r f belonging to the fixed Eulerian mesh feels a force from all surrounding membrane points. Further details of this numerical method applied to the simulation of vesicle dynamics (without viscosity contrast Λ = 1) can be found in our previous contribution (Kaoui et al. 2011).

Viscosity contrast
In the present paper, we consider the viscosities of the internal encapsulated fluid η int and of the external suspending fluid η ext to be different. We define the viscosity contrast as the ratio between these two viscosities: We are interested in cases where Λ = 1 like, for example, healthy red blood cells for which the hemoglobin solution (the internal fluid) is about five to seven times more viscous than the plasma (the external fluid). In 2D, a vesicle is presented by a discretized closed contour. It is then seen, from a geometry point of view, as a closed polygon. We find that a simple way to include the viscosity contrast in the algorithm is to attribute to each fluid node (r f ) either the viscosity, η int or η ext , depending on whether this point being located inside (Ω int ) or outside (Ω ext ) the polygon (i.e., the vesicle). To this purpose, we use the even-odd rule (O'Rourke 1998) to detect the location of each fluid node (of the Eulerian fluid mesh) at each time iteration. The even-odd rule consists of checking whether a point is located inside or outside a polygon by testing how many times a ray, starting from this point and going to any fixed direction, intersects the edges of the polygon. If the point in question is not on the boundary of the polygon, the number of the intersections is an even number if the point is outside, and it is odd if it is inside. This algorithm is also known as the crossing number algorithm, and it is largely used in computational geometry. Once the current location of a fluid node, with respect to the vesicle, is known, we assign to it the corresponding relaxation time For simplicity, τ ext is set to unity and τ int is calculated based on τ ext and Λ: τ int = Λ (τ ext − 1/2) + 1/2. This approach allows to describe the internal and external fluids using a single lattice Boltzmann equation. This requires fewer computational costs as compared to solving the internal and external fluids individually using two LB equations (Halliday et al. 2013). To account for the variation of the viscosity in space, we only change the relaxation time parameter depending on the location of the fluid nodes. This gives rise to a sharp jump of the viscosity across the membrane. This is generally more realistic than using a diffuse field for the viscosity as it is done in other front capturing methods such as phase-field or level-set methods (Biben and Misbah 2003;Beaucourt et al. 2004;Du et al. 2006;Maitre et al. 2012;Laadhari et al. 2012;Salac and Miksis 2012;Doyeux et al. 2013). In real world conditions, the membrane is a material barrier that separates two fluids with different viscosities. Our present method does not rely on diffuse scalar fields for the viscosity, and no additional field equation needs to be solved. Despite the simplicity of our method, it is capable of capturing the known physics of vesicles, as it is shown below. In Fig. 1, we summarize the main steps of our method.

A single particle
A single neutrally buoyant vesicle (the mass densities of its internal and external fluids are identical) is placed at mid-distance between the two parallel plates as shown in Fig. 2. A linear simple shear flow, with desired shear rate γ , is generated by moving the two plates in opposite directions. All results shown below have been obtained for a Reynolds number of 0 < Re = ργ R 2 0 /η ext < 1 and a capillary number of 0 < Ca = η ext γ R 3 0 /κ B < 10 , where R 0 is the effective vesicle radius. We set R 0 = 20 which allows to describe the deformation with good resolution and accuracy while still keeping the computational effort within acceptable limits. In 2D, R 0 = P /2π , where P is the vesicle perimeter. The swelling degree (sometimes called the reduced area) Δ = 4πA/P 2 (A is the vesicle area) is chosen as Δ = 0.8 in order to compare to the data in Beaucourt et al. (2004). The external viscosity η ext is held constant in all simulations. The parameter of interest in this paper is the viscosity contrast Λ, which is varied from 0.5 to 24 (0.5 < Λ < 24) by varying only the internal viscosity Fig. 1 The main steps of the numerical method η int . The size of the rectangular simulation box is taken as L x = 40 × R 0 (the length) and L y = 20 × R 0 (the height). L x is chosen to be large enough to minimize the effect of the periodic boundary conditions at the in-and outlet (an array of particles separated with longer distance L x to avoid interaction). In this way, we recover the limit of a single isolated vesicle.

Dynamical states under shear flow
A vesicle subjected to shear flow is known to undergo either tank-treading (steady motion) or tumbling (unsteady motion) (De Haas et al. 1997;Kantsler and Steinberg 2005;2006;Mader et al. 2006). These dynamical states are also observed for red blood cells (Fischer et al. 1978;Abkarian et al. 2007). In the tank-treading state (TT), the vesicle main axis assumes a steady inclination angle with respect to the flow direction while its membrane undergoes a tanktreading like motion. For the tumbling state (TB), the vesicle rotates as a solid elongated particle. The vesicle exhibits one of these dynamical states depending on three main dimensionless parameters Zhao and Shaqfeh 2011): the viscosity contrast Λ, the swelling degree Δ, and the capillary number Ca. Recently, we found that another fourth dimensionless parameter also controls the dynamics: the degree of confinement χ (the ratio of the vesicle diameter to the channel height) (Kaoui et al. 2012). χ should be small to avoid wall effects and therefore we use χ = 0.1 and 0.05.
We limit ourselves to the dynamical transition induced solely by varying the viscosity contrast Λ and reproduce numerically the tank-treading and the tumbling motions as depicted in the snapshots of Figs. 3 and 4. For smaller viscosity contrasts (e.g., Λ = 1), the vesicle tank-treads (Fig. 3). The shape of the vesicle and its orientation with respect to the flow direction are steady in time, while its membrane performs a tank-treading like motion: every point  on the membrane rotates around the center of mass of the vesicle. In Fig. 3, we see a marker point indicated on the membrane with a full black circle rotating in time around the vesicle center of mass. The tank-treading motion of the vesicle membrane generates a rotational flow in the internal encapsulated fluid. The applied shear flow reorients the vesicle and induces the motion of its membrane and of its internal fluid. The motion of the membrane on its turn disturbs the flow of the external fluid in its vicinity. By increasing Λ, and beyond a certain threshold, we trigger a dynamical transition from tank-treading to tumbling. The snapshots in Fig. 4 depict the tumbling motion of a vesicle with Λ = 15. The vesicle rotates in the same manner as a solid elongated particle. By increasing Λ, the internal medium of the vesicle becomes more and more viscous and therefore it dissipates more energy (Beaucourt et al. 2004).
In Fig. 5, we report the time evolution of the inclination angle θ (in degrees) for 1 ≤ Λ ≤ 13. This angle is computed as θ = arctan w y /w x , where w x and w y are the components of one of the eigenvectors of the tensor of inertia of the vesicle. In all simulations, θ is initially set to zero. This corresponds to the configuration of a vesicle being oriented parallel to the flow direction. Figure 5 exhibits all the quantitative features observed and known for vesicles when placed in a simple shear flow: for the tank-treading state, a transient time exists where θ increases until it reaches the steady state inclination angle θ * . This angle decreases with increasing the value of Λ. Above a critical value of Λ, the vesicle dynamics transits from tank-treading to tumbling. During the tumbling mode, θ varies in a periodic way similar to an arctan function of time. The period of the tumbling motion decreases with Λ. Vesicles with larger Λ rotate faster. This is a good benchmark that demonstrates the ability of our method to recover the TT and the TB states of a vesicle subjected to shear flow and also to trigger the transition from TT to TB by increasing Λ above a critical value, as it is known from experiments (Kantsler and Steinberg 2006;Mader et al. 2006).

Stability
After demonstrating the ability of our method to reproduce the known dynamics of a vesicle subjected to shear flow, we study its stability. As mentioned in "Vesicles (structure)" section, any numerical method used to simulate a vesicle has to fulfill two constraints (in 2D): the conservation of both the vesicle enclosed area A and of its perimeter P . A is expected to be constant in time since the internal enclosed fluid is considered to be an incompressible Newtonian fluid. However, due to numerical inaccuracies, we observe that A varies in time. In order to keep A constant, we include an additional membrane force proportional to the measured area variation and directed along the normal vector, κ A [A − A 0 ] n, where A is the actual area of the vesicle and A 0 its initial area to be taken as a reference value. κ A is an arbitrary numerical parameter chosen as to achieve a good conservation of A. It is taken to be as large as possible to maximize the impact of the correction force, but in the limit of still having numerically stable simulations. For the perimeter conservation, the tension σ introduced in "Vesicles (structure)" section is computed as σ (r, t) = κ S [Δs − Δs 0 ], where Δs is the actual distance (at time t) between the two adjacent membrane points and Δs 0 is its initial value taken as a referential value. κ S is another free numerical parameter which controls the way the vesicle perimeter is locally conserved. It is chosen to be as large as possible in order to achieve low variation but sufficiently low to keep the simulations numerically stable. In Fig. 6, we report the variation (in percent) of A and P in time associated to the data shown in Fig. 5. Both quantities are computed as the variation with respect to the initial values A 0 and P 0 of the area and the perimeter, respectively: (A − A 0 )/A 0 for the area and (P − P 0 )/P 0 for the perimeter. A good conservation of both quantities is a obtained with our present numerical method. For all used viscosity contrasts 1 ≤ Λ ≤ 13, the variation is less than 0.001 % for the area and less than 0.1 % for the perimeter. Figure 7 shows the evolution in time (scaled by the inverse of the shear rate 1/γ ) of the inclination angle θ of a tumbling vesicle with Λ = 15. The angle evolves in a periodic and stable way during many tumbling periods. This demonstrates the stability of the method, which is achieved despite the high value of Λ and the simplicity of the algorithm we used to include the viscosity contrast.

Comparison with other numerical methods
In Fig. 8, we report the steady inclination angle θ * of tank-treading vesicles versus their viscosity contrast Λ. All vesicles have the same swelling degree (Δ = 0.8). For the same flow conditions (Re = 0.05, Ca = 0.5) and only by increasing Λ, while holding all other parameters constant, we observe a decrease of θ * with Λ. Vesicles with more viscous internal fluid tend to align with the flow direction. θ * → 0 when Λ → Λ C , where Λ C is the critical viscosity contrast where the dynamical transition from tank-treading to tumbling takes place.
In the same figure, we report data from Ref. (Beaucourt et al. 2004) obtained by the boundary integral method (BIM)  (Beaucourt et al. 2004), the phase-field method (PFM) (Beaucourt et al. 2004) and the two-dimensional version of the Keller and Skalak Theory (KS2D) (Beaucourt et al. 2004;Keller and Skalak 1982). While there is a qualitative agreement, a quantitative agreement is not obtained as outlined in the text and the two-dimensional version of the Keller and Skalak theory (KS2D), respectively (Keller and Skalak 1982). Both sets of data (BIM and KS2D) are obtained under the assumption of an unbounded external flow (which corresponds to χ = 0) and in the exact Stokes limit (Re = 0). We notice that θ * obtained by our method (LBM) agrees qualitatively with BIM and KS2D since it decreases with Λ. However, it disagrees quantitatively. Our measured angles are larger than the ones obtained either by BIM or KS2D for the whole interval of 1 ≤ Λ ≤ 7. The difference becomes significant at larger Λ and the tank-treading-to-tumbling transition threshold is also shifted to higher values, i.e., from 3.7 (KS2D) and 4.5 (BIM) to 6.2 (LBM). The disagreement between data obtained with the three different methods can be explained as follows: (i) The KS2D model does not take into account the deformability of the particle, but assumes the vesicle to be a rigid non-deformable particle (this corresponds to Ca = 0). (ii) In BIM and LBM, the deformability and as such the capillary number is finite (Ca = 0.5). Second, in contrast to the BIM and KS2D (where only the Stokes equations are solved), the LBM recovers the solutions of the full Navier-Stokes equations, i.e., including the inertia terms: (∂u/∂t + u · ∇u). This brings the Reynolds number to the picture as an additional dimensionless parameter that also governs the vesicle dynamics. As such, the discrepancy observed in Fig. 8 between our data and the one in Ref. (Beaucourt et al. 2004) can to some extend be associated to the presence of inertia-despite the small value of Re we use. (iii) The KS2D and BIM data were obtained for an unbounded geometry, while our simulations use a channel with bounding walls at the top and bottom. Even by assuming a very weak confinement (χ = 0.05 or χ = 0.1), we expect the results to differ quantitatively from an unbounded flow (χ = 0).
The only available data for non-zero Reynolds number (Re = 10 −2 ) and in the presence of walls (χ = 0.1) is the phase-field (PFM) data presented in Beaucourt et al. (2004). Only a single data point is available and denoted by the upwards pointing triangle in Fig. 8. It is beyond the values obtained by other methods which do not account for bounding walls, but close to ours. This supports our argument that the walls are also responsible and an important factor to explain the difference between our data and the BIM. The impact of residual inertia and confinement is particularly important in 2D, where the hydrodynamic forces decay as ln(r)-in contrast to 3D simulations, where they decay at a much faster rate as 1/r (Kromkamp et al. 2006). Here, r is the distance separating a given point in the fluid from the point where a singular force is exerted. As such, increasing the size of the simulation box or decreasing Re further (by dropping the shear rate) is not sufficient to suppress the effect of residual inertia and confinement in 2D. For this reason, our code is more appropriate for simulations where confined geometries and inertia play a role. A systematic comparison between different methods in 3D is required in the future, but beyond the scope of the current paper. Following the above arguments, it is not surprising that our LBM simulations do not show a quantitative agreement with, e.g., BIM simulations. The effect of residual inertia, the presence of walls, and the choice of the capillary number strongly influence the dynamic behavior of viscous vesicles and as such can have a substantial impact on the tank-treading to tumbling transition.
Sensitivity to inertia Re = 0 Figure 9 depicts the dependence of θ * on Λ for different Reynolds numbers Re between 0.01 and 0.9. Re is varied by varying the shear rate γ , while keeping all other physical and numerical parameters unchanged. Re is taken between 0 and 1. Here, we do not study the inertia effects as in Refs. (Laadhari et al. 2012;Salac and Miksis 2012), but we investigate the contribution of the inertial terms present in our method. By varying γ alone, we consequently vary both parameters: Re and Ca. For Re < 0.1, the inclination angle is hardly sensitive to Re. However, for Re > 0.1, the θ * (Λ) curves shift upwards towards higher values of θ * when increasing Re. We also find a shift of the critical viscosity contrast Λ C needed for the tank-treading-to-tumbling transition towards larger values. A tumbling vesicle (at low Re) transits to tank-treading mode by increasing Re above a certain critical value. The same trend was also observed by Laadhari et al. (2012) and Salac and Miksis (2012) . 9 The steady inclination angle θ * as a function of the viscosity contrast Λ for tank-treading vesicles at different values of the Reynolds number Re. At low Re, the data is only weakly sensitive to inertia. At larger Re, the angle increases and the threshold of the tanktreading-to-tumbling transition shifts to larger values. Inertia tends to increase the steady inclination angle of tank-treading vesicles and to delay their transition to the tumbling mode. Here, Δ = 0.8 and χ = 0.1 tank-treading-to-tumbling transition changes its configuration with varying Re. For smaller Re, the average inclination angle during the tumbling motion is zero. At larger values, we observe that the average angle of the tumbling vesicles is not zero anymore, but adopts a finite positive value (e.g., 6 degrees for Re = 0.5).
In the two Refs. (Laadhari et al. 2012;Salac and Miksis 2012), the values of Λ and Re (e.g., Λ up to 50 and Re up to 100) are beyond what can be achieved in realistic laboratory experiments or in the microcirculatory system. For healthy red blood cells suspended in plasma, 7 < Λ < 13 and for lipid vesicles, Λ can be only increased experimentally up to 20 (Vitkova et al. 2009). If one considers giant unilamellar vesicles (with 10μm ≤ R 0 ≤ 100μm and 0.9 × 10 −19 ≤ κ B ≤ 20 × 10 19 J ) suspended in an aqueous solution (η ∼ 10 −3 kg/ms and ν ∼ 10 −6 m/s 2 ) and subjected to shear rates up to 100s −1 , the maximum Re that can be achieved is smaller than unity (Re < 1). This also explains the small range of Re chosen for this paper. Moreover, Re is related to Ca via the relationship which is independent of the flow and only depends on the properties of the suspending fluid and the vesicle. Smaller ratios Ca/Re are obtained for smaller and stiffer vesicles. Interestingly, in our simulation, by increasing Re, we increase Ca and so does θ * . This observed behavior contradicts the previous finding in ), where it was reported that θ * decreases with Ca for a tank-treading vesicle under shear flow. A finding that was confirmed also numerically later by Zhao and Shaqfeh (2011). However, this was obtained by the small deformation theory valid only in the Stokes limit (Re = 0) in contrast to the present work, where Re is not exactly zero. However, if we fix Re to a given value even not necessarily zero and we vary only Ca (by varying κ B ), we find again the same trend observed in Kaoui et al. (2009) and Zhao and Shaqfeh (2011): θ * decreases with Ca. The only way to vary Ca while keeping Re fixed is to vary the membrane rigidity κ B . If we vary solely the shear rate γ , we vary at the same time both Ca and Re because both parameters depend linearly on γ .

Multiple particles
So far we applied our numerical approach to the simple case of a single isolated vesicle placed in a simple shear flow.
Here, we present a more sophisticated application: multiple vesicles (with viscosity contrast) flowing in a microfluidic constriction (a complex geometry). We place six vesicles just before the inlet of the constriction. A flow generated by a body force advects the vesicles from the top to the bottom, see Figs. 10 and 11. The only possibility for the vesicles to pass through the constriction is to enter into it one after the other because the size of each vesicle is comparable to the width of the constriction (R 0 /W = 20/20 = 1). Once a vesicle enters and moves inside the constriction, its shape experiences large deformations (i.e., the vesicle is squeezed) due to the strong confinement and to the high flow rate. In Fig. 10, we show the case of non-viscous vesicles (Λ = 1). They deform and rearrange themselves in such a way that they enter and move in the constriction one after the other. During their motion in the constriction, they (where Λ = 1). A higher viscosity contrast affects the motion and the shape deformability of vesicles. It makes them stiffer and slower deform, and once they exit at the outlet, they tend to recover their equilibrium shape. This transition from a squeezed shape to the equilibrium shape is induced by the variation in the degree of confinement from higher R 0 /W = 1 (in the constriction) to lower R 0 /W = 0.36 (at the outlet). For a comparison and to qualitatively depict the effect of increasing the viscosity contrast, in Fig. 11, we show how viscous vesicles (with Λ = 15) behave. These vesicles are less deformable. First, we notice that it takes some time for the first vesicle to enter the constriction and once it has entered it moves slowly. At the outlet, more time is required to recover the equilibrium shape as compared to the nonviscous vesicles (Fig. 10). Such a microfluidic geometry is used for medical diagnosis of diseases that affect the RBCs deformability (Shelby et al. 2003). For example, it is known that malaria makes RBCs less deformable than healthy ones. The time required for the cells to recover their equilibrium shape at the outlet of the constriction can be used as a smart tool to evaluate the stiffness of diseased RBCs.

Conclusions
We introduced a numerical method to simulate in 2D the dynamics of vesicles with viscosity contrast (their internal and external fluids have different dynamics viscosities). The flow of the internal and the external fluids is computed using the lattice Boltzmann method, and the fluidvesicle (fluid-structure) two-way coupling is accomplished by the immersed boundary method. The viscosity contrast is included in our numerical method by attributing different viscosities to a given fluid node depending on its location with respect to the vesicle. We use a simple algorithm inspired from the field of computational geometry which allows to test whether a fluid node is inside or outside the vesicle. We demonstrated the capabilities of our approach to reproduce the known dynamics of a vesicle with viscosity contrast subjected to shear flow. The dynamical transition from tank-treading (a liquid-like motion) to tumbling (a solid-like motion) is triggered by increasing the viscosity of the internal fluid beyond a threshold. We examined the stability of our code, and we compared our simulation data against theory and data obtained by other numerical methods. The effect of inertia was discussed briefly. Even if the Reynolds number is taken smaller than unity, it does not mean the inertia effects are absent. With our code, we capture similar inertia-dependent dynamical behavior as reported also numerically by others (Salac and Miksis 2012;Laadhari et al. 2012;Luo et al. 2013). However, it has to be pointed out that these recent numerical works are not yet accompanied by experiments confirming the numerical predictions. Our suggested method can handle also multiple viscous vesicles flowing in complex geometries. As a proof of concept, we performed simulations of six viscous vesicles flowing through a microfluidic constriction.
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.