Numerical Study of a Confined Vesicle in Shear Flow at Finite Temperature

Abstract

The dynamics and rheology of a vesicle confined in a channel under shear flow are studied at finite temperature. The effect of finite temperature on vesicle motion and system viscosity is investigated. A two-dimensional numerical model, which includes thermal fluctuations and is based on a combination of molecular dynamics and mesoscopic hydrodynamics, is used to perform a detailed analysis in a wide range of the Peclet numbers (the ratio of the shear rate to the rotational diffusion coefficient). The suspension viscosity is found to be a monotonous increasing function of the viscosity contrast (the ratio of the viscosity of the encapsulated fluid to that of the surrounding fluid) both in the tank-treading and the tumbling regime due to the interplay of different temperature-depending mechanisms. Thermal effects induce shape and inclination fluctuations of the vesicle which also experiences Brownian diffusion across the channel increasing the viscosity. These effects reduce when increasing the Peclet number.

1. Introduction

Suspensions of soft particles such as droplets, vesicles, and capsules are ubiquitous in relevant applications in biology, medicine, and engineering. Studying their dynamics in flow is challenging, since shapes are not fixed, as in the case of rigid objects, but depend dynamically on the interplay between fluid stresses and interfacial forces. The interfacial forces are directly related to the nature of the considered particles: The surface tension for droplets, the membrane bending rigidity for vesicles, and additionally the membrane shear elasticity for capsules. This calls for separate investigations of the various systems.

Vesicles are small volumes of fluid embedded in a lipid bi-layer membrane, in solution with either the same or different fluid. The dynamical and rheological properties of their suspensions in flow have attracted a lot of theoretical and experimental interest, as comprehensively reviewed in Refs. [1,2,3,4]. A consensus has been reached concerning the dynamical regimes in shear flow. In dilute solution, vesicles can show tank-treading (TT), tumbling (TU), and vacillating-breathing (VB) (also called trembling or swinging) motion, depending on the shear rate and the viscosity contrast $\lambda ={\eta }_{out}/{\eta }_{in}$, where ${\eta }_{in}$ and ${\eta }_{out}$ are the viscosities of the inner and outer fluids, respectively. TT and TU occur at low and high $\lambda$, respectively, while VB appears for strong flows when vesicle deformation affects its dynamics [5,6,7,8,9,10,11,12,13,14,15].

On the other hand, the rheology of single vesicle suspensions is still a matter of debate. Indeed, different behaviors of the intrinsic viscosity ${\eta }_I=(\eta -{\eta }_{out})/\left({\eta }_{out}\varphi \right)$, where $\eta$ is the effective system viscosity and $\varphi$ the vesicle concentration, as a function of the viscosity contrast have been observed in experimental, theoretical, and numerical studies. In the case of very dilute suspensions of quasi-spherical vesicles, it was shown analytically [16,17] that the intrinsic viscosity decreases with the viscosity contrast $\lambda$ in the TT regime, reaching a minimum at the TT-to-TU transition, and then grows with $\lambda$ in the TU regime. Experimental investigations do not provide conclusive results. A good agreement with the theoretical prediction was found in Ref. [18], while an increase in ${\eta }_I$ with $\lambda$ for $\lambda <1$ was observed in Ref. [19]. These discrepancies might be due to the difficulty in preparing monodisperse suspensions as well as to the fact that viscosity measurements require volume fractions $\varphi \sim$ 5–10%, thus making the extrapolation to the dilute limit difficult [19]. Numerical models differ mainly in the lack or presence of thermal noise. In the former case, it was found in two-dimensional models that the intrinsic viscosity follows the theoretical prediction both in the very dilute [20,21,22] and in the dilute case [21,23,24]. A similar dependence on the viscosity contrast was found also in a three-dimensional model [25]. The only available numerical model with thermal fluctuations [26] shows that ${\eta }_I$ is an increasing function of $\lambda$, in agreement with the experiments of Ref. [19].

The numerical model of Ref. [26], which comprises both thermal membrane undulations and thermal noise [26], is adopted here to perform a detailed study of a confined vesicle in shear flow at a finite temperature. The results of this model yielded very good agreement with the experimental results in describing the collision process of two vesicles [19] and the flow field of a single vesicle in shear flow [27]. The system is studied in two dimensions at a fixed shear rate in a wide range of the Peclet number $Pe$—the ratio of the shear rate to the rotational diffusion coefficient— differently from other theoretical and numerical studies where $Pe=\infty$. We aim at elucidating the role played by thermal fluctuations in influencing both the vesicle dynamics and, consequently, the system viscosity in the TT and TU regimes. The reason for considering a very dilute solution is twofold. On one hand, this is in line with the hypothesis of an extremely dilute suspension used in the theoretical model [16,17], and, on the other hand, hydrodynamic and steric interactions between vesicles can be ruled out.

The paper is organized as follows. Section 2 presents the numerical model. Results are illustrated in Section 3. A detailed discussion of our findings about the effects of thermal noise is presented in Section 4, including a comparison with previous studies. Finally, conclusions are presented in Section 5.

2. The Model

A two-dimensional fluid made of $N_s$ point-like particles of mass m is considered. The particle positions ${\mathbf{r}}_i\left(t\right)$ and velocities ${\mathbf{v}}_i\left(t\right)$, $i=1,2,\dots ,N_s$, at time t are continuous variables. We employ the multi-particle collisions (MPC) dynamics approach, in which time evolution occurs via iterative propagations and collisions [28,29,30,31]. In the first streaming step, particles are ballistically streamed for a time interval $\Delta t_s$

(1)${\mathbf{r}}_i(t+\Delta t_s)={\mathbf{r}}_i\left(t\right)+{\mathbf{v}}_i\left(t\right)\Delta t_{s}i=1,\dots ,N_s.$

In the subsequent collision step, the system is divided into square cells of mesh size a where an instantaneous multi-particle collision occurs, which changes particle velocities as

(2)${\mathbf{v}}_i^{new}={\mathbf{v}}_c^G+{\mathbf{v}}_i^{ran}-\sum _{j\in cell}{\mathbf{v}}_j^{ran}/N_c+{\mathsf{\Pi }}^{-1}\sum _{j\in cell}m\left[{\mathbf{r}}_{j,c}\times ({\mathbf{v}}_j-{\mathbf{v}}_j^{ran})\right]\times {\mathbf{r}}_{i,c}i=1,\dots ,N_s$

where ${\mathbf{v}}_c^G$ is the center-of-mass velocity of all particles in the cell, ${\mathbf{v}}_i^{ran}$ is a velocity taken from a Maxwell–Boltzmann distribution, $N_c$ is the number of particles in the cell, $\mathsf{\Pi }$ and ${\mathbf{r}}_{i,c}$ are the moment-of-inertia tensor and the position relative to the center of mass of the particles in the cell, respectively. This dynamic conserves both local linear and angular momentum [32,33] and keeps the temperature constant [34]. The viscosity of the fluid is given by [35]

(3)$\eta =\frac{m}{\Delta t_s}\left[{\left(\frac{l}{a}\right)}^2\left(\frac{n^2}{n-1}-\frac{n}{2}\right)+\frac{1}{24}\left(n-\frac{7}{5}\right)\right]$

n being the average number of particles per cell, $l=\Delta t_s\sqrt{k_{B}T/m}$ the mean-free path, and $k_{B}T$ the thermal energy. The system of size $L_x\times L_y$ is confined between two horizontal walls sliding along the x direction with velocities $v_{wall}$ and $-v_{wall}$. Periodic boundary conditions (BC) are used along the x direction. Bounce-back BC are enforced at walls [36] obtaining a linear flow profile $(u_x,u_y)=(\dot{\gamma }y,0)$ with shear rate $\dot{\gamma }=2v_{wall}/L_y$.

The vesicle membrane is modeled as a chain of $N_p$ beads of mass $m_p$ connected to form a closed ring with an average bond length $r_0$. Neighboring beads interact via a harmonic potential

(4)$U_{bond}={\kappa }_h\sum _{i=1}^{N_p}\frac{\left(\right|{\mathbf{r}}_i-{\mathbf{r}}_{i-1}{|-r_0)}^2}{2r_0^2}$

where ${\kappa }_h$ is the spring constant and ${\mathbf{r}}_i$ is the position vector of the i-th bead. This ensures the conservation of the membrane length. Shapes and fluctuations are controlled by the bending potential

(5)$U_{bend}=\frac{\kappa }{r_0}\sum _{i=1}^{N_p}(1-cos{\beta }_i),$

where $\kappa$ is the bending rigidity and ${\beta }_i$ is the angle between two consecutive bonds. Finally, the internal area S is kept close to the target area $S_0$ of the vesicle by using a quadratic constraint-potential with a compression modulus ${\kappa }_S$ [26]

(6)$U_{area}={\kappa }_S\frac{{(S-S_0)}^2}{2r_0^4}.$

Newton’s equations of motions of beads are integrated by using the velocity-Verlet algorithm with time step $\Delta t_p$ [37].

In order to describe the coupling of solvent particles with the vesicle, each bead is treated as a “rough" hard disk radius $r_v$ [26,38,39]. The value of $r_v$ is set so that disks overlap and achieve full coverage of the membrane. Scattering takes place when a solvent particle i and a disk j overlap while moving towards each other so that both the conditions $|{\mathbf{r}}_j-{\mathbf{r}}_i|<r_v$ and $({\mathbf{r}}_j-{\mathbf{r}}_i)·({\mathbf{v}}_j-{\mathbf{v}}_i)<0$ are fulfilled. A second disk $k=j\pm 1$, connected to the j-th one and characterized by the smallest distance from the solvent particle i, is then selected. The angular velocity

(7)$\mathsf{\Omega }={\mathsf{\Pi }}^{-1}\sum _{l=i,j,k}{m}_l{\mathbf{r}}_{l,c}\times {\mathbf{v}}_l$

and the center of mass velocity ${\mathbf{v}}^G$ of the $i,j,k$-particle system are computed, with ${\mathbf{r}}_{l,c}$ being the position relative to the center of mass. The updated values of the velocities are given by

(8)${\mathbf{v}}_l^{new}=2({\mathbf{v}}^G+\mathsf{\Omega }\times {\mathbf{r}}_{l,c})-{\mathbf{v}}_{l}l=i,j,k$

which guarantees linear and angular momenta conservation [14]. The collision step (2) is then performed for those fluid particles which did not interact with the membrane in order to avoid multiple collisions with the same membrane disk in the following iterations. Disks interact with lateral walls also by implementing bounce-back scattering. The numerical implementation of the algorithm is outlined in Appendix A.

Inertial effects, which are experimentally irrelevant due to the small flow velocities, are made negligible in the simulations by making the Reynolds number $Re=\dot{\gamma }\rho R_0^2/{\eta }_{out}$, with mass density $\rho$, which is very small. Other relevant dimensionless quantities are the reduced area $S^{*}=S_0/\pi R_0^2$, where $R_0=L_0/2\pi$ is the vesicle radius with $L_0$ the vesicle contour length, and the reduced shear rate ${\dot{\gamma }}^{*}=\dot{\gamma }{\tau }_c$, where ${\tau }_c={\eta }_{out}{R}_0^3/\kappa$ is the relaxation time of the vesicle. The viscosity contrast can be approximated as $\lambda \simeq m_{in}/m_{out}$ within the present model [32] (the subscripts $out/in$ refer to quantities outside/inside the vesicle). We use the following $L_x=18.95R_0$, $L_y=5.79R_0$ with $R_0=7.6a$. Finally, we set $m_{in}$ so as to obtain $0.1\le \lambda \le 15.0$, $m_p=3m_{out}$, $\Delta t_s/\Delta t_p=64$, $N_p=480$, $r_v=r_0=a/10$, ${\kappa }_S=4\times {10}^{-4}{k}_{B}T$, ${\kappa }_h=3\times {10}^2{k}_{B}T$. The setting of parameters is implemented to obtain $Re<0.15$, the Mach number $Ma=v_{wall}/c_s<0.25$, where $c_s=\sqrt{2k_{B}T/m_{out}}$ is the speed of sound, to reduce compressibility effects [40], and ${\dot{\gamma }}^{*}=1.0$ in all the cases. The value of the reduced shear rate ${\dot{\gamma }}^{*}$ is comparable to those used in other studies [20,21,23,25] and provides access to the TT and TU regimes by varying the viscosity contrast.

The importance of thermal fluctuations depends on the the rotational Peclet number $Pe=\dot{\gamma }/D_r$. The rotational diffusion coefficient $D_r$ is given by $D_r=k_{B}T/\zeta$ and by employing the rotational friction coefficient $\zeta$ of a circle, the Peclet number can be written as $Pe=4\pi {\dot{\gamma }}^{*}\kappa /\left(k_{B}T{R}_0\right)$. In the following, the Peclet number will be changed by considering the values $\kappa /\left(k_{B}T{R}_0\right)=6.58,65.8,164.5$, corresponding to $Pe=82,821,2041$, respectively, while keeping fixed the value of ${\dot{\gamma }}^{*}$. The present study focuses on the dynamics and rheology of a sheared vesicle at finite values of $Pe$. Indeed, in previous studies of Refs. [20,21,23,24,25] it was assumed that $Pe=\infty$, thus neglecting the role of thermal fluctuations.

3. Results

We consider very dilute suspensions with a single vesicle for two values of the reduced area $S^{*}=0.80,0.95$ corresponding effectively to volume fractions $\varphi =0.023,0.028$, respectively.

In Figure 1, the instantaneous intrinsic viscosity ${\eta }_I$ is shown as a function of time for different values of viscosity contrast $\lambda$, bending rigidity $\kappa$, and reduced area $S^{*}$. The viscosity $\eta$ is computed as $\eta ={\sigma }_{xy}/\dot{\gamma }$ where ${\sigma }_{xy}$ is the $xy$ component of the stress tensor at walls [41]. In the MPC model, the stress ${\sigma }_{xy}$ has a contribution in the streaming step, ${\sigma }_{xy}^s$, proportional to the flux of the x-momentum crossing the walls, and a second contribution in the collision step, ${\sigma }_{xy}^c$, due to the multi-particle collision with virtual wall particles (see Appendix A). In two-dimensional simulations the streaming contribution is [42]

(9)${\sigma }_{xy}^s=\frac{m}{L_x\Delta t_s}\sum _{i=1}^{N_s}[v_{x,i}^{\prime }\left(t_b\right)-v_{x,i}\left(t_b\right)],$

where $t_b$ ($t\le t_b\le t+\Delta t_s$) is the time when particle i bounces back from the wall, $v_{x,i}^{\prime }\left(t_b\right)$ and $v_{x,i}\left(t_b\right)$ are the velocities just after and before the collision with the wall, respectively, and $N_s$ is the number of particles hitting one of the walls in the time interval $[t,t+\Delta t_s]$. The collision contribution is [42]

(10)${\sigma }_{xy}^c=\frac{m}{L_x\Delta t_s}\sum _{i=1}^{N_c}[v_{x,i}^{\prime }(t+\Delta t_s)-v_{x,i}(t+\Delta t_s)],$

where $N_c$ is the number of particles with multi-particle collision with virtual wall particles, while $v_{x,i}^{\prime }(t+\Delta t_s)$ and $v_{x,i}(t+\Delta t_s)$ are the velocities of particle i after and before the collision step, respectively.

After a transient period, when the vesicle moves from the initial position towards the center of the channel attaining its steady state, ${\eta }_I$ fluctuates around average values up to the longest simulated times, which are more than two orders of magnitude larger than the vesicle relaxation time ${\tau }_c$.

The values $⟨{\eta }_I⟩$ of the intrinsic viscosity, time-averaged in the steady state, are reported in Figure 2 as a function of $\lambda$. It appears that $⟨{\eta }_I⟩$ is an increasing function of $\lambda$ for the used values of the reduced area, bending energy, and temperature, in agreement with our previous results [26,39]. In the Keller–Skalak theory [5], where thermal fluctuations are ignored, the sharp TT-to-TU transition occurs at ${\lambda }_c\simeq 3.7$ for $S^{*}=0.80$ and at ${\lambda }_c\simeq 6.5$ for $S^{*}=0.95$. However, finite temperature broadens the TT-to-TU transition [14]. In the TU regime at higher values of $\lambda$, the growth of $⟨{\eta }_I⟩$ is steeper. A decrease in the intrinsic viscosity in the TT regime followed by its growth in the TU regime, as theoretically predicted in Refs. [16,17] and observed in simulations without thermal fluctuations [20,21,23,24,25], is not found in our model. The effect of increasing the bending energy is to reduce the value of the intrinsic viscosity without changing the monotonic dependence on the viscosity contrast. This effect seems to be triggered by the Peclet number as will be discussed later.

In order to clarify the observed behavior of $⟨{\eta }_I⟩$, the vesicle dynamics were investigated in more detail by monitoring the temporal evolution of several quantities. The inclination angle $\mathsf{\Theta }$, describing the angle between the x direction and the long main axis of the vesicle, can be used to discriminate between the TT and the TU states. In the former case, $\mathsf{\Theta }$ reaches a steady value, while in the latter case, $\mathsf{\Theta }$ varies periodically in time. In Figure 3, the inclination angle is shown as a function of time. For low values of $\lambda$ the vesicle performs tank-treading motion and the inclination angle fluctuates around a steady value. In contrast, without thermal fluctuations [21,23] the inclination angle is constant in the TT regime after the initial transient. When increasing the viscosity contrast, some tumbling events appear, which become predominant for the highest value of $\lambda$.

The time-averaged values $⟨\mathsf{\Theta }⟩$ are depicted in Figure 4, together with the root-mean-square (rms) fluctuation values ${\sigma }_{\mathsf{\Theta }}=\sqrt{⟨{(\Delta \mathsf{\Theta })}^2⟩}$. The transition from the TT to the TU regime, which is characterized by going from values $⟨\mathsf{\Theta }⟩>0$ to $⟨\mathsf{\Theta }⟩\simeq 0$, is broader for the smallest values of the bending rigidity, and becomes sharper when increasing the ratio $\kappa /\left(k_{B}T{R}_0\right)$. The fluctuations ${\sigma }_{\mathsf{\Theta }}$ reduce with $Pe$ in the TT regime, as theoretically predicted [38], and show an opposite trend with increasing viscosity contrast.

From the gyration tensor of the vesicle, the two eigenvalues ${\mathsf{\Lambda }}_M$ and ${\mathsf{\Lambda }}_m$ with ${\mathsf{\Lambda }}_M>{\mathsf{\Lambda }}_m$ are extracted and the asphericity $A={[({\mathsf{\Lambda }}_M-{\mathsf{\Lambda }}_m)/({\mathsf{\Lambda }}_M+{\mathsf{\Lambda }}_m)]}^2$ is computed.

The values of A as a function of time are shown in Figure 5 and the time-averages $⟨A⟩$ as a function of the viscosity contrast in Figure 6.

$⟨A⟩$ is constant in the TT regime and decreases when approaching the TU regime, showing that the vesicle becomes more rounded when the inner fluid is more viscous. Additionally, $⟨A⟩$ is smaller for the lower value of bending rigidity and does not change significantly going from TT to TU regime for the highest value of the bending rigidity. In the case of the quasi-circular vesicle a non-monotonic behavior of $⟨A⟩$ with the bending rigidity can be observed in the TT regime. The average values $⟨\sqrt{{\mathsf{\Lambda }}_M}⟩$ and $⟨\sqrt{{\mathsf{\Lambda }}_m}⟩$, which provide an estimate of the vesicle semi-axes, are plotted in Figure 7 as a function of the viscosity contrast to demonstrate how the vesicle becomes more rounded when increasing $\lambda$ for fixed ${\dot{\gamma }}^{*}$.

It can be seen that $⟨\sqrt{{\mathsf{\Lambda }}_M}⟩$ decreases and $⟨\sqrt{{\mathsf{\Lambda }}_m}⟩$ increases as functions of $\lambda$. The relative change of the average eigenvalues, from the TT to the TU regime, is larger at $\kappa /\left(k_{B}T{R}_0\right)=65.8$, while it is negligible for the highest value of the bending rigidity. The rms fluctuation values ${\sigma }_M=\sqrt{⟨{(\Delta \sqrt{{\mathsf{\Lambda }}_M})}^2⟩}$ and ${\sigma }_m=\sqrt{⟨{(\Delta \sqrt{{\mathsf{\Lambda }}_m})}^2⟩}$ are reported in Figure 8 as functions of $\lambda$. In all the cases, the values of the rms fluctuations are constant in the TT regime and increase when entering the TU regime. Moreover, ${\sigma }_M$ and ${\sigma }_m$ decrease when increasing the Peclet number.

The time behavior of the vertical position $y_{cm}$ of the vesicle center of mass displays Brownian diffusion across the channel width up to the longest simulated time, as can be seen in Figure 9.

The vesicle does not span the whole channel cross-section due to the lift force which pushes it far from the walls [14]. In previous studies [21,23], where thermal noise is absent, vesicles move along the center line of the channel without lateral displacement and with a regular arrangement in the TT steady state, in two or three files at higher concentrations [43,44]. It was later found that there is a critical viscosity contrast above which the vesicle can be either placed along the center line or off-centered without lateral wandering [45]. The rms fluctuation values ${\sigma }_{cm}=\sqrt{⟨{(\Delta y_{cm})}^2⟩}$ are reported in Figure 10. For the lowest values of the bending rigidity it is evident that ${\sigma }_{cm}$ increases with the viscosity ratio $\lambda$ due to the more circular shape, while this trend is less pronounced for further increases in $\kappa /\left(k_{B}T{R}_0\right)$.

Moreover, a reduction in the values of ${\sigma }_{cm}$ can be observed when increasing the bending rigidity with no significant dependence on the reduced area $S^{*}$. In the TT regime it becomes ${\sigma }_{cm}/R_0\simeq \sqrt{\left(k_{B}T{R}_0\right)/\kappa }\propto \sqrt{1/Pe}$ for the explored range of bending rigidities. The term $\sqrt{\left(k_{B}T{R}_0\right)/\kappa }$ is the rms value of the vesicle deformation amplitude [38].

Finally, the average configurations of the vesicle are presented in Figure 11 for reduced area $S^{*}=0.80,0.95$, bending rigidity $\kappa /\left(k_{B}T{R}_0\right)=6.58,65.8$, and two values of the viscosity contrast. The shapes are obtained by averaging in time and space, in the vesicle eigenvector reference frame, the positions of membrane beads in circular sectors of width $\pi /45$ radians. This visualizes how the vesicle becomes more rounded from the TT to the TU regime in the case with $\kappa /\left(k_{B}T{R}_0\right)=65.8$ at $S^{*}=0.95$. The reduction in the asphericity is less appreciable in the other cases.

4. Discussion

We can now relate the observed behavior of the intrinsic viscosity $⟨{\eta }_I⟩$ in Figure 2 to the changes in vesicle shape and diffusion. We think that the monotonic growth of $⟨{\eta }_I⟩$ is due to the interplay of several mechanisms. As previously observed in Ref. [19], shape fluctuations favor energy dissipation that increases $⟨{\eta }_I⟩$, while alignment with the flow direction causes a decrease in $⟨{\eta }_I⟩$ with increasing viscosity contrast. The vesicle becomes more rounded with increasing $\lambda$ as indicated by the average asphericity. As a consequence, the vesicle experiences a larger resistance to the flow with the tilt angle approaching the values $\pi /4$. This counteracts the reduction due to the decrease in the average inclination angle when approaching the TT-to-TU transition. The most relevant effect due to thermal noise of the fluid is that the vesicle is not located at the center of the channel, but wanders across it due to fluctuation-induced Brownian diffusion (the possible influence of this effect on the intrinsic viscosity was already mentioned in Ref. [23]). This implies that the vesicle can never move along the centerline of the channel, which is the state of minimum dissipation when thermal effects are neglected [23]. The amplitude of this lateral motion is quantified by ${\sigma }_{cm}$, which grows with an increasing viscosity ratio for the lowest value of the bending rigidity. Since the vesicle becomes closer to the walls, a larger resistance of the vesicle to the flow might be induced, similarly to what happens for colloids whose effective diffusion coefficient reduces close to a wall [46]. This effect would contribute to the increase in the $⟨{\eta }_I⟩$ even in the TT regime. We remark that since it results to be ${\sigma }_{cm}\propto \sqrt{1/Pe}$, as previously found, much higher values of the Peclet number are required in order to access a regime where ${\sigma }_{cm}\simeq 0$ to ignore thermal fluctuations.

The outlined picture persists when increasing the bending rigidity $\kappa$, when the value of $⟨{\eta }_I⟩$ is reduced but its $\lambda$-dependence is not affected. Similar values of $⟨{\eta }_I⟩$ are observed for the highest bending rigidity where the TT-to-TU transition is sharper and the vesicle becomes more rigid, as observed in the values of the average asphericity and of the rms fluctuations ${\sigma }_M$ and ${\sigma }_m$ which hardly change with $\lambda$. In the TT regime, the effect of increasing the bending rigidity is to reduce the average inclination angle $⟨\mathsf{\Theta }⟩$, its variance ${\sigma }_{\mathsf{\Theta }}$, and ${\sigma }_{cm}$ with respect to the case with the lowest bending rigidity, while the vesicle appears to be less circular. As a consequence the vesicle has less resistance to the flow, which explains the reduction in $⟨{\eta }_I⟩$ when compared to lower values of $\kappa$. In the TU regime, the difference in the average asphericity for the three values of the bending rigidity diminishes, causing a reduction in the difference of the average intrinsic viscosities.

To complete our discussion, we note that it was argued in Ref. [23] that the monotonic behavior of $⟨{\eta }_I⟩$ might be due to measurements performed in short transient regimes; however, as here shown, this is not the case. Moreover, our results do not depend on the choice either of the channel length $L_x/R_0=19$ or of the degree of confinement $2R_0/L_y=0.35$, as suggested in Refs. [21,23]. Indeed, these two values are intermediate between the ones used in those studies [21,23] where the non-monotonic behavior of the intrinsic viscosity was observed without thermal fluctuations.

5. Conclusions

We believe that the monotonic growth of $⟨{\eta }_I⟩$ has to be related to the presence of thermal fluctuations that are missing in other models. This effect persists up to the highest Peclet number of about $2\times {10}^3$. In a simplified stochastic three-dimensional model of vesicles in shear flow [7], it was shown that thermal fluctuations cannot be neglected at up to $Pe=1.2\times {10}^3$. Much higher values of $Pe$ are required, as previously discussed, in order to ignore thermal fluctuations. Finally, we add that the relevance of thermal noise in the vesicle dynamics was demonstrated also for the VB regime in numerical [11,14], theoretical [47,48], and experimental studies [9,49].

Data Availability Statement

Data are available upon reasonable request.

Acknowledgments

A.L. wishes to thank G. Gompper for useful discussions and hospitality at Forschungszentrum Jülich. This work was performed under the auspices of GNFM-INdAM.

Conflicts of Interest

The author declares no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Time behavior of the intrinsic viscosity [Forumla omitted. See PDF.] (values are averaged over time intervals of duration [Forumla omitted. See PDF.] to smooth out noise) at [Forumla omitted. See PDF.] (upper row) and [Forumla omitted. See PDF.] (lower row) for [Forumla omitted. See PDF.] (left), [Forumla omitted. See PDF.] (right) with [Forumla omitted. See PDF.] 1 (full green line), 7 (dashed blue line), 15 (dot-dashed red line).

Figure 2. Average values of the intrinsic viscosity [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.] (right) for [Forumla omitted. See PDF.] (•), 65.8 (▲), 164.5 (⋆). Full lines are guides to the eye. The tank-treading to tumbling transition occurs at [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] and at [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] in the Keller–Skalak theory [5] and is marked by the dashed vertical lines. Error bars are given by the root-mean-square fluctuation values of the intrinsic viscosity.

Figure 3. Time behavior of the inclination angle [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.] (right) for [Forumla omitted. See PDF.] with [Forumla omitted. See PDF.] 1 (full green line), 7 (dashed blue line), 15 (dot-dashed red line).

Figure 4. Average values of the inclination angle [Forumla omitted. See PDF.] (filled symbols) and its rms fluctuation values [Forumla omitted. See PDF.] (empty symbols) as a function of [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.] (right) for [Forumla omitted. See PDF.] (•), 65.8 (▲), 164.5 (⋆). The tank-treading to tumbling transition in the Keller–Skalak theory [5] is marked by the dashed vertical lines.

Figure 5. Time behavior of the asphericity A (values are sampled every [Forumla omitted. See PDF.]) at [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.] (right) for [Forumla omitted. See PDF.] with [Forumla omitted. See PDF.] 1 (full green line), 7 (dashed blue line), 15 (dot-dashed red line).

Figure 6. Average values of the asphericity [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.] (right) for [Forumla omitted. See PDF.] (•), 65.8 (▲), 164.5 (⋆). The tank-treading to tumbling transition in the Keller–Skalak theory [5] is marked by the dashed vertical lines. Error bars are given by the root-mean-square fluctuation values of the asphericity.

Figure 7. Average values of the gyration tensor eigenvalues [Forumla omitted. See PDF.] (filled symbols) and [Forumla omitted. See PDF.] (empty symbols) as a function of [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] (□), [Forumla omitted. See PDF.] (+), and at [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] (•), 65.8 (▲), 164.5 (⋆). The tank-treading to tumbling transition in the Keller–Skalak theory [5] is marked by the dashed ([Forumla omitted. See PDF.]) and full ([Forumla omitted. See PDF.]) vertical lines.

Figure 8. Rms fluctuation values [Forumla omitted. See PDF.] (filled symbols) and [Forumla omitted. See PDF.] (empty symbols) of the the gyration tensor eigenvalues of Figure 7 as a function of [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.] (right) for [Forumla omitted. See PDF.] (•), 65.8 (▲), 164.5 (⋆). The tank-treading to tumbling transition in the Keller–Skalak theory [5] is marked by the dashed vertical lines.

Figure 9. Time behavior of the vertical position [Forumla omitted. See PDF.] of the vesicle center of mass (values are sampled every [Forumla omitted. See PDF.]) at [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.] (right) for [Forumla omitted. See PDF.] with [Forumla omitted. See PDF.] 1 (full green line), 7 (dashed blue line), 15 (dot-dashed red line). The horizontal full line denotes the center of the channel.

Figure 10. Rms fluctuation values [Forumla omitted. See PDF.] of the vertical position of the vesicle center of mass of Figure 9 as a function of [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.] (right) for [Forumla omitted. See PDF.] (•), 65.8 (▲), 164.5 (⋆). The tank-treading to tumbling transition in the Keller–Skalak theory [5] is marked by the dashed vertical lines.

Figure 11. Vesicle average configurations at [Forumla omitted. See PDF.] (upper row) and [Forumla omitted. See PDF.] (lower row) for [Forumla omitted. See PDF.] (left), [Forumla omitted. See PDF.] (right) with [Forumla omitted. See PDF.]. The full line represents the unit circle as a reference.

Appendix A

Here, the numerical implementation of the algorithm is described for a system bounded by two moving solid walls. In the presence of walls, the system consists of fluid particles and virtual particles. The fluid particles represent the solvent, while the virtual particles are required to impose no-slip conditions at the walls. First, the positions and velocities, for both the solvent particles and the vesicle beads, are initialized. The fluid particles are distributed uniformly inside the system with average number n of particles per cell. An extra layer of collision cells is required next to the walls to enforce boundary conditions. The virtual particles are also uniformly distributed with the same number density. All the fluid real particles and beads are initialized with velocities sampled from the Maxwell–Boltzmann distribution with variances $k_{B}T/m$ and $k_{B}T/m_p$, respectively, and zero mean. The velocities of virtual particles are from the Maxwell-–Boltzmann distribution with variances $k_{B}T/m$ and average $\pm \frac{1}{2}\dot{\gamma }{L}_y$. The initial linear and angular momentum are removed from each cell and from all the beads, and the velocities are rescaled to set the temperature to the value T.

At each time step Newton’s equations of motion for beads are integrated by means of the velocity-Verlet algorithm with time step $\Delta t_p$ [37]. Every $\Delta t_s/\Delta t_p$ time steps the MPC algorithm and the solvent-vesicle collisions are performed in the following way:

• All the solvent particles are streamed according to Equation (1). Particles crossing walls undergo bounce-back collisions changing their velocities as ${\mathbf{v}}_i\to 2{\mathbf{v}}_{wall}-{\mathbf{v}}_i$ where ${\mathbf{v}}_{wall}$ and $-{\mathbf{v}}_{wall}$ are the wall velocities with ${\mathbf{v}}_{wall}=(v_{wall},0)$.

• The solvent particles and the beads which overlap, are looked for and their velocities are modified according to Equation (8).

• Galilean invariance is violated when the mean-free path l is much smaller than the cell size a. To restore the Galilean invariance [50], all the fluid particles are moved by a random vector $\mathbf{s}$ as ${\mathbf{r}}_i\to {\mathbf{r}}_i+\mathbf{s}$. The components of this random vector are drawn from a uniform distribution in the interval $[-a/2,a/2]$.

• All solvent particles are sorted in respective cells and cell-level quantities are calculated.

• The velocities ${\mathbf{v}}_i$ of fluid particles not scattering with the vesicle, are updated according to Equation (2). The virtual particles are assigned a new random velocity.

• All fluid particles are shifted back to their original position as ${\mathbf{r}}_i\to {\mathbf{r}}_i-\mathbf{s}$.

References

1. Vlahovska, P.M.; Podgorski, T.; Misbah, C. Vesicles and red blood cells: From individual dynamics to rheology. C. R. Phys.; 2009; 10, 775. [DOI: https://dx.doi.org/10.1016/j.crhy.2009.10.001]

2. Abreu, D.; Levant, M.; Steinberg, V.; Seifert, U. Fluid vesicles in flow. Adv. Colloid Interface Sci.; 2014; 208, 129. [DOI: https://dx.doi.org/10.1016/j.cis.2014.02.004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24630339]

3. Winkler, R.G.; Fedosov, D.A.; Gompper, G. Dynamical and rheological properties of soft colloid suspensions. Curr. Opin. Colloid Interface Sci.; 2014; 19, 594.

4. Barthès-Biesel, D. Motion and deformation of elastic capsules and vesicles in flow. Annu. Rev. Fluid Mech.; 2016; 48, 25. [DOI: https://dx.doi.org/10.1146/annurev-fluid-122414-034345]

5. Keller, S.R.; Skalak, R. Motion of a tank-treading ellipsoidal particle in a shear flow. J. Fluid. Mech.; 1982; 120, 27. [DOI: https://dx.doi.org/10.1017/S0022112082002651]

6. Noguchi, H.; Gompper, G. Fluid vesicles with viscous membranes in shear flow. Phys. Rev. Lett.; 2004; 93, 258102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.93.258102]

7. Noguchi, H.; Gompper, G. Dynamics of fluid vesicles in shear flow: Effect of membrane viscosity and thermal fluctuations. Phys. Rev. E; 2005; 72, 011901. [DOI: https://dx.doi.org/10.1103/PhysRevE.72.011901]

8. Kantsler, V.; Steinberg, V. Orientation and dynamics of a vesicle in tank-treading motion in shear flow. Phys. Rev. Lett.; 2005; 95, 258101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.95.258101]

9. Kantsler, V.; Steinberg, V. Transition to tumbling and two regimes of tumbling motion of a vesicle in shear flow. Phys. Rev. Lett.; 2006; 96, 036001. [DOI: https://dx.doi.org/10.1103/PhysRevLett.96.036001]

10. Misbah, C. Vacillating breathing and tumbling of vesicles under shear flow. Phys. Rev. Lett.; 2006; 96, 028104. [DOI: https://dx.doi.org/10.1103/PhysRevLett.96.028104]

11. Noguchi, H.; Gompper, G. Swinging and tumbling of fluid vesicles in shear flow. Phys. Rev. Lett.; 2007; 98, 128103. [DOI: https://dx.doi.org/10.1103/PhysRevLett.98.128103] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17501159]

12. Lebedev, V.V.; Turitsyn, K.S.; Vergeles, S.S. Dynamics of nearly spherical vesicles in an external flow. Phys. Rev. Lett.; 2007; 99, 218101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.99.218101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18233260]

13. Vlahovska, P.M.; Gracia, R.S. Dynamics of a viscous vesicle in linear flows. Phys. Rev. E; 2007; 75, 016313. [DOI: https://dx.doi.org/10.1103/PhysRevE.75.016313] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17358259]

14. Messlinger, S.; Schmidt, B.; Noguchi, H.; Gompper, G. Dynamical regimes and hydrodynamic lift of viscous vesicles under shear. Phys. Rev. E; 2009; 80, 011901. [DOI: https://dx.doi.org/10.1103/PhysRevE.80.011901]

15. Zhao, H.; Shaqfeh, E.S.G. The dynamics of a vesicle in simple shear flow. J. Fluid Mech.; 2011; 674, 578. [DOI: https://dx.doi.org/10.1017/S0022112011000115]

16. Danker, G.; Misbah, C. Rheology of a dilute suspension of vesicles. Phys. Rev. Lett.; 2007; 98, 088104. [DOI: https://dx.doi.org/10.1103/PhysRevLett.98.088104]

17. Danker, G.; Biben, T.; Podgorski, T.; Verdier, C.; Misbah, C. Dynamics and rheology of a dilute suspension of vesicles: Higher-order theory. Phys. Rev. E; 2007; 76, 041905. [DOI: https://dx.doi.org/10.1103/PhysRevE.76.041905]

18. Vitkova, V.; Mader, M.A.; Polack, B.; Misbah, C.; Podgorski, T. Micro-macro link in rheology of erythrocyte and vesicle suspensions. Biophys. J.; 2008; 95, L33. [DOI: https://dx.doi.org/10.1529/biophysj.108.138826]

19. Kantsler, V.; Segre, E.; Steinberg, V. Dynamics of interacting vesicles and rheology of vesicle suspension in shear flow. EPL; 2008; 82, 58005. [DOI: https://dx.doi.org/10.1209/0295-5075/82/58005]

20. Ghigliotti, G.; Biben, T.; Misbah, C. Rheology of a dilute two-dimensional suspension of vesicles. J. Fluid Mech.; 2010; 653, 489. [DOI: https://dx.doi.org/10.1017/S0022112010000431]

21. Kaoui, B.; Jonk, R.J.W.; Harting, J. Interplay between microdynamics and macrorheology in vesicle suspensions. Soft Matter; 2014; 10, 4735. [DOI: https://dx.doi.org/10.1039/C4SM00563E] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24851823]

22. Nait-Ouhra, A.; Farutin, A.; Ez-Zahraouy, H.; Benyoussef, A.; Misbah, C. Rheology of a confined vesicle suspension. Phys. Rev. Fluids; 2019; 4, 103602. [DOI: https://dx.doi.org/10.1103/PhysRevFluids.4.103602]

23. Thiébaud, M.; Misbah, C. Rheology of a vesicle suspension with finite concentration: A numerical study. Phys. Rev. E; 2013; 88, 062707. [DOI: https://dx.doi.org/10.1103/PhysRevE.88.062707] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24483486]

24. Rahimian, A.; Veerapaneni, S.K.; Biros, G. Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method. J. Comput. Phys.; 2010; 229, 6466. [DOI: https://dx.doi.org/10.1016/j.jcp.2010.05.006]

25. Zhao, H.; Shaqfeh, E. The dynamics of a non-dilute vesicle suspension in a simple shear flow. J. Fluid Mech.; 2013; 725, 709. [DOI: https://dx.doi.org/10.1017/jfm.2013.207]

26. Lamura, A.; Gompper, G. Dynamics and rheology of vesicle suspensions in wall-bounded shear flow. EPL; 2013; 102, 28004. [DOI: https://dx.doi.org/10.1209/0295-5075/102/28004]

27. Afik, A.; Lamura, A.; Steinberg, V. Long-range hydrodynamic effect due to a single vesicle in linear flow. EPL; 2016; 113, 38003. [DOI: https://dx.doi.org/10.1209/0295-5075/113/38003]

28. Malevanets, A.; Kapral, R. Mesoscopic model for solvent dynamics. J. Chem. Phys.; 1999; 110, 8605. [DOI: https://dx.doi.org/10.1063/1.478857]

29. Malevanets, A.; Kapral, R. Solute molecular dynamics in a mesoscale solvent. J. Chem. Phys.; 2000; 112, 7260. [DOI: https://dx.doi.org/10.1063/1.481289]

30. Kapral, R. Multiparticle Collision Dynamics: Simulation of Complex Systems on Mesoscales. Adv. Chem. Phys.; 2008; 140, 89.

31. Gompper, G.; Ihle, T.; Kroll, D.M.; Winkler, R.G. Multi-Particle Collision Dynamics: A Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids. Adv. Polym. Sci.; 2009; 221, 1.

32. Noguchi, H.; Kikuchi, N.; Gompper, G. Particle-based mesoscale hydrodynamic techniques. Europhys. Lett.; 2007; 78, 10005. [DOI: https://dx.doi.org/10.1209/0295-5075/78/10005]

33. Götze, I.O.; Noguchi, H.; Gompper, G. Relevance of angular momentum conservation in mesoscale hydrodynamics simulations. Phys. Rev. E; 2007; 76, 046705. [DOI: https://dx.doi.org/10.1103/PhysRevE.76.046705] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17995137]

34. Allahyarov, A.; Gompper, G. Mesoscopic solvent simulations: Multiparticle-collision dynamics of three-dimensional flows. Phys. Rev. E; 2002; 66, 036702. [DOI: https://dx.doi.org/10.1103/PhysRevE.66.036702] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12366293]

35. Noguchi, H.; Gompper, G. Transport coefficients of off-lattice mesoscale-hydrodynamics simulation techniques. Phys. Rev. E; 2008; 78, 016706. [DOI: https://dx.doi.org/10.1103/PhysRevE.78.016706]

36. Lamura, A.; Gompper, G.; Ihle, T.; Kroll, D.M. Multi-particle collision dynamics: Flow around a circular and a square cylinder. Europhys. Lett.; 2001; 56, 319. [DOI: https://dx.doi.org/10.1209/epl/i2001-00522-9]

37. Allen, M.P.; Tildesley, D.J. Computer Simulation of Liquids; Clarendon Press: Oxford, UK, 1987.

38. Finken, R.; Lamura, A.; Seifert, U.; Gompper, G. Two-dimensional fluctuating vesicles in linear shear flow. Eur. Phys. J. E; 2008; 25, 309. [DOI: https://dx.doi.org/10.1140/epje/i2007-10299-7]

39. Lamura, A.; Gompper, G. Rheological properties of sheared vesicle and cell suspensions. Procedia IUTAM; 2015; 16, 3. [DOI: https://dx.doi.org/10.1016/j.piutam.2015.03.002]

40. Lamura, A.; Gompper, G. Numerical study of the flow around a cylinder using multi-particle collision dynamics. Eur. Phys. J. E; 2002; 9, 477. [DOI: https://dx.doi.org/10.1140/epje/i2002-10107-0]

41. Mewis, J.; Wagner, N.J. Colloidal Suspension Rheology; Cambridge University Press: Cambridge, UK, 2012.

42. Tao, Y.-G.; Götze, I.O.; Gompper, G. Multiparticle collision dynamics modeling of viscoelastic fluids. J. Chem. Phys.; 2008; 128, 144902. [DOI: https://dx.doi.org/10.1063/1.2850082]

43. Thiébaud, M.; Shen, Z.; Harting, J.; Misbah, C. Prediction of anomalous blood viscosity in confined shear flow. Phys. Rev. Lett.; 2014; 112, 238304. [DOI: https://dx.doi.org/10.1103/PhysRevLett.112.238304] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24972235]

44. Shen, Z.; Farutin, A.; Thiébaud, M.; Misbah, C. Interaction and rheology of vesicle suspensions in confined shear flow. Phys. Rev. Fluids; 2017; 2, 103101. [DOI: https://dx.doi.org/10.1103/PhysRevFluids.2.103101]

45. Nait-Ouhra, A.; Guckenberger, A.; Farutin, A.; Ez-Zahraouy, H.; Benyoussef, A.; Gekle, S.; Misbah, C. Lateral vesicle migration in a bounded shear flow: Viscosity contrast leads to off-centered solutions. Phys. Rev. Fluids; 2018; 3, 123601. [DOI: https://dx.doi.org/10.1103/PhysRevFluids.3.123601]

46. Brenner, H. The slow motion of a sphere through a viscous fluid towards a plane surface. Chem. Eng. Sci.; 1961; 16, 242. [DOI: https://dx.doi.org/10.1016/0009-2509(61)80035-3]

47. Abreu, D.; Seifert, U. Effect of thermal noise on vesicles and capsules in shear flow. Phys. Rev. E; 2012; 86, 010902. [DOI: https://dx.doi.org/10.1103/PhysRevE.86.010902] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23005361]

48. Abreu, D.; Seifert, U. Noisy nonlinear dynamics of vesicles in flow. Phys. Rev. Lett.; 2013; 110, 238103. [DOI: https://dx.doi.org/10.1103/PhysRevLett.110.238103] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25167533]

49. Levant, M.; Steinberg, V. Amplification of thermal noise by vesicle dynamics. Phys. Rev. Lett.; 2012; 109, 268103. [DOI: https://dx.doi.org/10.1103/PhysRevLett.109.268103]

50. Ihle, T.; Kroll, D.M. Stochastic rotation dynamics: A Galilean-invariant mesoscopic model for fluid flow. Phys. Rev. E; 2001; 63, 020201(R). [DOI: https://dx.doi.org/10.1103/PhysRevE.63.020201]