2. Red blood cell modelling strategies

We have implemented the vesicle and capsule models in our platform, which have been previously described in Boedec et al. (Reference Boedec, Leonetti and Jaeger2017) and Lyu et al. (Reference Lyu, Chen, Boedec, Leonetti and Jaeger2021) for use in unbounded and confined spaces, respectively. The vesicle model accounts for the fluid nature of the membrane and enforces surface incompressibility by projecting the three-dimensional (3-D) velocity field onto a space with zero surface divergence (Boedec, Leonetti & Jaeger Reference Boedec, Leonetti and Jaeger2011). In addition, the isogeoparametric representation of the geometry (Cottrell, Hughes & Bazilevs Reference Cottrell, Hughes and Bazilevs2009), which has at least $C^1$ (i.e. continuous first derivatives) regularity, allows for the incorporation of the Helfrich bending energy using a rigorous weak mathematical formulation. Meanwhile, the capsule model can utilise several models of polymerised membrane behaviour and incorporate bending elasticity through the Helfrich formulation developed for the vesicle model or other common forms of thin-shell modelling. Although the vesicle and capsule models can incorporate shear and dilatation surface viscosity developed for a surfactant drop (Gounley et al. Reference Gounley, Boedec, Jaeger and Leonetti2016), we will not be utilising this option for our objective of this paper, which is to assess the purely elastic aspects of RBC mechanics. As the details of these models are available elsewhere, we will focus solely on their key elements in this section.

Our isogeometric representation is obtained using the Loop subdivision method (Loop Reference Loop1987; Cirak, Ortiz & Schröder Reference Cirak, Ortiz and Schröder2000; Cirak & Ortiz Reference Cirak and Ortiz2001). To generate the control mesh for a vesicle or capsule, we start with an icosahedron and perform $N$ Loop subdivisions. The value $N=0$ corresponds to the initial icosahedron, which has $N_{elem}=20$ equilateral triangular elements that correspond to 20 Loop elements. With each refinement level, each element is subdivided into four equivalent equilateral triangles, and a node is inserted at the centre of each edge. The spatial position of all nodes in the refined mesh is determined using Loop's rules, which ensure zero geometric representation error. The level $N$ generates a mesh with $N_{elem}=20 \times 4^N$. We typically achieve satisfactory accuracy using $N=3$, resulting in $N_{elem}=1280$. Note that the twelve nodes of the $N=0$ level have fifth-order connectivity, while all other nodes introduced afterwards have sixth-order connectivity. These twelve nodes are called singular, and they are an inevitable trace of the icosahedron. The regularity of the Loop approximation is only $C^1$ at these nodes, while it is $C^2$ everywhere else. It is worth highlighting that our isogeometric finite-element method, which relies on the Loop subdivision, diverges from traditional shell modelling techniques commonly used in structural mechanics. The advantage of our approach is its ability to circumvent the shear-locking phenomenon typically observed in conventional methods (Cirak et al. Reference Cirak, Ortiz and Schröder2000; Cirak & Ortiz Reference Cirak and Ortiz2001).

The boundary element method (Pozrikidis Reference Pozrikidis1992, Reference Pozrikidis2002) is used with the surface mesh described earlier to obtain the velocity of the lipid medium at each node. We can then determine the velocity at any point of an element using Loop's interpolation functions. The normal component of this velocity field determines the evolution of the surface position between two time steps in the simulation. Meanwhile, the movement of the nodes tangentially to the membrane is identified with the kinematics of a cytoskeleton. This tangential movement is determined by equilibrating the elastic forces exerted by the cytoskeleton with the viscous friction forces of the points of attachment of the cytoskeleton in the lipid bilayer. The viscous friction force is proportional to the velocity difference between the bilayer and the cytoskeleton. For the unconfined case and without considering a viscosity contrast to simplify the presentation, the system can be formulated as follows (Lyu et al. Reference Lyu, Chen, Boedec, Leonetti and Jaeger2018):

(2.1)\begin{gather} \boldsymbol{v}^{BL} = \boldsymbol{v}^{3D} = \boldsymbol{\mathsf{P}} \left( \boldsymbol{v}^{\infty} + \boldsymbol{\mathsf{G}} \boldsymbol{f}^{RBC} \right) \quad \mbox{with} \ \boldsymbol{f}^{RBC} = \boldsymbol{f}^{BL} + \boldsymbol{f}^{SC}, \end{gather}

(2.2)\begin{gather}\boldsymbol{x}(t+\mathrm{d}t) = \boldsymbol{x}(t) + \boldsymbol{v}^{SC}\,\mathrm{d}t \quad \mbox{with} \ \boldsymbol{v}^{SC} = \boldsymbol{v}^{BL} + \boldsymbol{v}^{SC/BL}. \end{gather}

Position and velocity fields of the RBC membrane are denoted by $\boldsymbol {x}$ and $\boldsymbol {v}$, respectively. Superscripts BL and SC are used to distinguish between the lipid bilayer and the spectrin cytoskeleton, respectively, while $3D$ is for the bulk flow and $\infty$ for the imposed background one. The velocity difference between the cytoskeleton and the lipid bilayer $\boldsymbol {v}^{SC/BL} = \boldsymbol {v}^{SC} - \boldsymbol {v}^{BL}$ is a tangential vector. The Green's operator associated with the Stokeslet is represented by $\boldsymbol{\mathsf{G}}$, and the projection operator on a subspace with zero surface divergence is denoted by $\boldsymbol{\mathsf{P}}$. This operator ensures that the surface tension $\gamma$, which is the Lagrange multiplier of the surface incompressibility constraint, satisfies the condition that the surface divergence of the velocity field $\boldsymbol {v}^{BL}$ is zero. The surface densities of force induced by the lipid bilayer and the cytoskeleton correspond to $\boldsymbol {f}^{BL}$ and $\boldsymbol {f}^{SC}$, respectively, and their sum $\boldsymbol {f}^{RBC}$ corresponds to the surface force density exerted by the entire RBC membrane on the ambient fluids. The matrix expression of the system (2.1)–(2.2) and the principle of the solution algorithms are presented in detail in Boedec et al. (Reference Boedec, Leonetti and Jaeger2017) and Lyu et al. (Reference Lyu, Chen, Boedec, Leonetti and Jaeger2018, Reference Lyu, Chen, Boedec, Leonetti and Jaeger2021). While we encourage interested readers to refer to these sources, we note that the rest of the paper is accessible without prior knowledge of them.

Using this validated approach, which has been successfully applied in our previous studies of vesicles and capsules, we have developed several strategies for modelling RBCs. These strategies encompass single-layer vesicle and capsule models, as well as double-layer vesicle–capsule and capsule–capsule models. To aid in visual comprehension, figure 1 presents a schematic depiction of the surface densities relevant to the RBC models. In the subsequent sections, we provide a detailed outline of the implementation for each of these models. To aid in the comparison of the results across all the studies, the same colour coding is utilised in all figures, where the vesicle is represented in black, the capsule in blue, the vesicle–capsule in red and the capsule–capsule in green.

Figure 1.

Schematic representation of the RBC models illustrating the relevant surface force densities on the respective layers. (a) Single-layer vesicle model with out-of-plane elasticity ($\boldsymbol {f}^{e\perp }$) and surface incompressibility ($\boldsymbol {f}^{\gamma }$). (b) Single-layer capsule model with in-plane ($\boldsymbol {f}^{e\parallel }$) and out-of-plane elasticities (surface incompressibility is approximatively considered). (c) Double-layer vesicle–capsule model with in-plane and out-of-plane elasticities, along with surface incompressibility. (d) Double-layer capsule–capsule model with in-plane and out-of-plane elasticities (surface incompressibility is approximatively considered). All models employ a single mesh, with the distinction that, in double-layer models, $\boldsymbol {v}^{BL} \neq \boldsymbol {v}^{SC}$ due to the differing assumed tangential kinematics of the two layers.

Throughout this paper, we use dimensionless variables, denoted by a star symbol. Lengths are expressed in units of a volume-based radius $R\equiv [3V/(4{\rm \pi} )]^{1/3}$, where $V\simeq 94\ \mathrm {\mu }{\rm m}^3$ represents the enclosed volume of RBCs. The surface area $A$ is approximately $135\ \mathrm {\mu }{\rm m}^2$, giving a reduced volume $v\equiv 3\sqrt {4{\rm \pi} }VA^{-3/2} = 0.64$ (Evans & Skalak Reference Evans and Skalak1980). Time is scaled by $\eta _{ext}R/\mu _s$, where $\eta _{ext}$ is the viscosity of the suspending fluid and $\mu _s$ is the surface shear modulus. A summary of the mechanical properties utilised in the simulations is provided in table 1. These data are based on averages of recognised values for a healthy RBC (see, e.g. table 1 in Levant & Steinberg Reference Levant and Steinberg2016). The table also displays the distribution of properties between the cytoskeleton and the lipid bilayer in the modelling strategies that distinguish them.

Table 1.

Mechanical properties.

The reference shape, whether quasi-spherical or biconcave, represents the shape relative to which in-plane and out-of-plane deformations are defined. In simpler terms, when the cytoskeleton (Skalak capsule) adopts the reference shape, the in-plane deformations are zero. It is important to clarify that all our simulations begin with an RBC already in the biconcave shape (i.e. $v=0.64$), which we obtain directly from an analytical expression (Evans & Skalak Reference Evans and Skalak1980) and consistently used in our prior work, as outlined in Lyu et al. (Reference Lyu, Chen, Boedec, Leonetti and Jaeger2018). The reference shapes we consider include the biconcave shape itself, a quasi-spherical shape ($v=0.96$) with zero spontaneous curvature and a quasi-spherical shape ($v=0.96$) with positive curvature ($C_0^* = C_0/R = 4$) (Peng et al. Reference Peng, Mashayekh and Zhu2014). If not specified, we have used quasi-spherical with zero spontaneous curvature. The quasi-spherical reference shape, achieved through a relaxation process, corresponds to the deflated shape of a vesicle with a reduced volume of $v=0.96$.

2.1. Single-layer vesicle strategy (colour code $=$ black)

The vesicle model does not consider a separate contribution of the cytoskeleton (i.e. $\boldsymbol {f}^{SC} = \boldsymbol {0}$), while

(2.3)\begin{gather} \boldsymbol{f}^{BL}(\boldsymbol{x}) = \boldsymbol{f}^{\gamma}(\boldsymbol{x}) + \boldsymbol{f}^{e\perp}(\boldsymbol{x}) ={-}\frac{\delta \mathcal{W}^{BL}}{\delta \boldsymbol{x}} ={-}\frac{\delta ( \mathcal{W}^{\gamma} + \mathcal{W}^{H} )}{\delta \boldsymbol{x}}, \end{gather}

(2.4a–c)\begin{gather}\mathcal{W}^{\gamma} = \int_S \gamma \,\mathrm{d} S, \quad \mathcal{W}^{H} = \int_S w^H \,\mathrm{d}S, \quad w^H = \frac{k_b}{2} (2 H + C_0 )^2. \end{gather}

The $\boldsymbol {f}^{e\perp }$ contribution represents the out-of-plane elastic forces, namely the bending forces induced by the Helfrich surface energy density $w^H$, which depend on the mean curvature $H$ and the spontaneous curvature $C_0$. The spontaneous curvature $C_0$ is defined as $C_0 = - 2H_0$, where $H_0$ is the mean curvature of the spontaneous or reference shape of the lipid bilayer. For zero spontaneous curvature, the Helfrich energy reduces to the simple expression $w^H = 2 k_b H^2$, where $k_b=2.4\times 10^{-19}$ J is the bending modulus. We omit the contribution of the Gaussian curvature, which can occur in $w^H$ but is not useful due to the Gauss–Bonnet theorem, as we are not considering a change in surface topology. The contribution of $\boldsymbol {f}^{\gamma }$ corresponds to the forces induced by the surface tension $\gamma$. For a constant tension, it results in an out-of-plane force contribution, but it can also produce in-plane forces if the tension varies, as $\boldsymbol {f}^{\gamma } = 2 \gamma H \boldsymbol {n} + \boldsymbol {\nabla }_s \gamma$, where $\boldsymbol {n}$ is the normal vector to the surface and $\boldsymbol {\nabla }_s= (\boldsymbol{\mathsf{I}} - \boldsymbol {n} \boldsymbol {n}) \boldsymbol {\cdot } \boldsymbol {\nabla }$ is the surface gradient operator, with $\boldsymbol{\mathsf{I}}$ being the identity in $\mathbb {R}^3$.

We want to emphasise that in the vesicle model, surface incompressibility, expressed as $\boldsymbol {\nabla }_s \boldsymbol {\cdot }\boldsymbol {v}^{BL} = 0$, is enforced through a projection algorithm, specifically, the projection operator $\boldsymbol{\mathsf{P}}$ as defined in (2.1). This projection introduces the surface tension parameter $\gamma$, acting as a Lagrange multiplier for the surface incompressibility constraint. While this method effectively ensures both local and global surface area conservation, it does come at a significant computational cost, as previously discussed (Boedec et al. Reference Boedec, Leonetti and Jaeger2017).

2.2. Single-layer capsule strategy (colour code $=$ )

The capsule model does not consider a separate contribution of the cytoskeleton (i.e. $\boldsymbol {f}^{SC} = \boldsymbol {0}$). While bending is still modelled using the Helfrich energy, the in-plane elasticity of the cytoskeleton is nicely represented by a capsule. On the other hand, fluidity is lost and surface incompressibility can only be approximated since the projector $\boldsymbol{\mathsf{P}}$ is replaced by the identity. Thus, we obtain the capsule model with

(2.5)\begin{equation} \boldsymbol{f}^{BL}(\boldsymbol{x}) = \boldsymbol{f}^{e \parallel}(\boldsymbol{x}) + \boldsymbol{f}^{e\perp}(\boldsymbol{x}) ={-}\frac{\delta \mathcal{W}^{BL}}{\delta \boldsymbol{x}} ={-}\frac{\delta ( \mathcal{W}^{SK} + \mathcal{W}^{H} )}{\delta \boldsymbol{x}} , \end{equation}

where $\boldsymbol {f}^{e\parallel }$ represents the in-plane elastic forces deduced from a polymerised membrane strain energy, which is usually defined on the reference configuration ($S^0$). Surface deformations are computed relative to this reference configuration (Boedec et al. Reference Boedec, Leonetti and Jaeger2017, (24)). To model the strain-hardening behaviour of the RBC membrane, we typically use the expression proposed by Skalak et al. (Reference Skalak, Tozeren, Zarda and Chien1973), which is given by

(2.6a,b)\begin{equation} \mathcal{W}^{SK} = \int_{S^0} w^{SK} \,\mathrm{d} S^0, \quad w^{SK} = \frac{G_s}{4} ( I_1^2 + 2I_1 - 2I_2 + C^{SK}I_2^2 ), \end{equation}

where the two invariants $I_1$ and $I_2$ are expressed as functions of the two main strains $\lambda _1$ and $\lambda _2$. The resistance to shear is controlled by the coefficient $G_s$, which is the surface shear modulus $\mu _s = 6\ \mathrm {\mu }{\rm N}\ {\rm m}^{-1}$, and the resistance to local area variation is controlled by the coefficient $C^{SK}$, which is the ratio of the area dilatation modulus to the shear modulus. The limiting case of surface incompressibility is obtained by making $C^{SK}$ tend towards infinity, where $I_2$ cancels. In practice, however, the numerical problem becomes too steep when $C^{SK}$ exceeds a hundred (Dodson & Dimitrakopoulos Reference Dodson and Dimitrakopoulos2010). For long-term simulations, we add a global area conservation constraint to prevent a strong drift of the RBC surface area. This constraint results in a contribution of the same type as $\mathcal {W}^{\gamma }$, with the surface tension $\gamma$ being replaced by a constant tension $k_s (S - S^0)/S^0$. In all our simulations, we used $C^{SK} = 80$ in combination with $k_s = 10^{3}\ \mathrm {\mu }{\rm N}\ {\rm m}^{-1}$ (Sigüenza, Mendez & Nicoud Reference Sigüenza, Mendez and Nicoud2017) when surface incompressibility has to be considered. Our extensive investigation has revealed that the contribution of $k_s$ remains negligible, as local area preservation is consistently maintained throughout our study.

2.3. Double-layer vesicle–capsule strategy (colour code $=$ )

The vesicle–capsule model is obtained by combining the expression (2.3) of $\boldsymbol {f}^{BL}$ and

(2.7)\begin{equation} \boldsymbol{f}^{SC}(\boldsymbol{x}) = \boldsymbol{f}^{e \parallel}(\boldsymbol{x}) ={-}\frac{\delta \mathcal{W}^{SC}}{\delta \boldsymbol{x}} ={-} \frac{\delta \mathcal{W}^{SK}}{\delta \boldsymbol{x}}. \end{equation}

In contrast to the vesicle model, the contribution from the cytoskeleton $\boldsymbol {f}^{SC}$ is not set to zero. The incompressibility constraint is still rigorously ensured by the $\boldsymbol{\mathsf{P}}$ projector. The coefficient $C^{SK}$ in Skalak's law is then fixed to represent the contribution of the cytoskeleton only, which is negligible compared with the resistance to expansion provided by the lipid bilayer. The dilatation modulus of the cytoskeleton was determined experimentally to be within the range of $1\unicode{x2013}10\ \mathrm {\mu }{\rm N}\ {\rm m}^{-1}$ (Lenormand et al. Reference Lenormand, Hénon, Richert, Siméon and Gallet2001, Reference Lenormand, Hénon, Richert, Siméon and Gallet2003; Kim et al. Reference Kim, Kim and Park2012), with a theoretical value for a regular hexagonal lattice being twice the shear modulus (Discher, Mohandas & Evans Reference Discher, Mohandas and Evans1994; Hansen et al. Reference Hansen, Skalak, Chien and Hoger1996). Thus, we set the expansion rigidity of the capsule representing the cytoskeleton to $C^{SK} = 2$ in our double-layer strategy simulations. In other words, the cytoskeleton only provides resistance to shear deformation, with $G_s = 6\ \mathrm {\mu }{\rm N}\ {\rm m}^{-1}$ and $C^{SK}=2$.

2.5. Membrane elastic force determination

The boundary conditions between the two-layer structure are implicit. Because of the single mesh representation, the two structures are constrained to have the same motion in the normal direction. In the tangent plane, sliding is permitted and limited exclusively by tangential friction forces. Therefore, the double-layer models depend on the frictional coupling between the lipid bilayer and cytoskeleton through junction protein complexes. Let $\boldsymbol {f}^{BL/SC}$ be the surface density of frictional force exerted by the lipid bilayer on the cytoskeleton. By the action–reaction principle, $\boldsymbol {f}^{SC/BL} = -\boldsymbol {f}^{BL/SC}$. The static equilibrium of the two structures can be expressed as

(2.8a,b)\begin{equation} \boldsymbol{f}^{SC} + \boldsymbol{f}^{BL/SC} = 0, \quad \boldsymbol{f}^{BL} + \boldsymbol{f}^{SC/BL} + \boldsymbol{f}^{ext} = 0, \end{equation}

where $\boldsymbol {f}^{ext}$ represents the action of the ambient fluids on the RBC membrane, which is in equilibrium with the forces induced by the deformation of the RBC, namely

(2.9)\begin{equation} \boldsymbol{f}^{ext} ={-}\boldsymbol{f}^{RBC} ={-}(\boldsymbol{f}^{BL} + \boldsymbol{f}^{SC}), \end{equation}

consistent with (2.1).

The finite-element method is used to determine the external forces through a weak formulation using the virtual work principle, for a given membrane configuration

(2.10)\begin{gather} \delta \mathcal{W}^{ext} + \delta \mathcal{W}^{int} = 0, \end{gather}

(2.11)\begin{gather} \delta \mathcal{W}^{ext}= \int_S \boldsymbol{f}^{ext} \boldsymbol{\cdot} \delta\boldsymbol{x}\,\mathrm{d}S ={-} \int_S \boldsymbol{f}^{RBC} \boldsymbol{\cdot} \delta\boldsymbol{x}\,\mathrm{d}S , \end{gather}

(2.12)\begin{align} \delta \mathcal{W}^{int}&={-}( \delta \mathcal{W}^{BL} + \delta \mathcal{W}^{SC} ) = \int_S -\frac{\delta \mathcal{W}^{BL}}{\delta \boldsymbol{x}} \boldsymbol{\cdot} \delta\boldsymbol{x}\,\mathrm{d}S + \int_S -\frac{\delta \mathcal{W}^{SC}}{\delta \boldsymbol{x}} \boldsymbol{\cdot} \delta\boldsymbol{x}\,\mathrm{d}S \nonumber\\ & = \int_S \boldsymbol{f}^{BL} \boldsymbol{\cdot} \delta\boldsymbol{x}\,\mathrm{d}S + \int_S \boldsymbol{f}^{SC} \boldsymbol{\cdot} \delta\boldsymbol{x}\,\mathrm{d}S. \end{align}

The internal virtual work $\delta \mathcal {W}^{int}$ depends on the surface stress and strain tensors, which can be determined based on the membrane configuration, see Boedec et al. (Reference Boedec, Leonetti and Jaeger2017) for details

(2.13)\begin{equation} \delta \mathcal{W}^{int} ={-}\int_S \left[ \sigma^{\alpha\beta} \delta(E_{\alpha\beta}) + \mu^{\alpha\beta} \delta(B_{\alpha\beta})\right]\mathrm{d}S = \int_S \left[\frac{1}{2} \sigma^{\alpha\beta} \delta(a_{\alpha\beta}) + \mu^{\alpha\beta} \delta(b_{\alpha\beta})\right]\mathrm{d}S. \end{equation}

Here, $\sigma _{\alpha \beta }$ and $\mu _{\alpha \beta }$ correspond to the stresses induced by membrane and bending strains, respectively. The tensor of components $E_{\alpha \beta }$ with $2E_{\alpha \beta } = a_{\alpha \beta }-a^{0}_{\alpha \beta }$ is the surface Green–Lagrange strain tensor and the tensor of components $B_{\alpha \beta } = b_{\alpha \beta }-b^{0}_{\alpha \beta }$ is the bending equivalent, where superscript 0 refers to the reference configuration. The $a_{\alpha \beta }$ are components of the metric tensor $\boldsymbol{\mathsf{a}}$ that correspond to the identity operator in the tangent plane $\boldsymbol{\mathsf{I}}^S = \boldsymbol{\mathsf{I}} - \boldsymbol {n} \boldsymbol {n}$. Its determinant $a = \mathrm {det}(\boldsymbol{\mathsf{a}})$ expresses the surface element as a function of the surface parameterisation $(s^1, s^2)$ as $\mathrm {d}S = \sqrt {a} \,\mathrm {d}s^1 \,\mathrm {d}s^2$.

2.6. Cytoskeleton drag force determination

We define $C_{f_{JC}}$ as the mean friction coefficient of a protein junction complex in the lipid bilayer. The frictional force generated by the movement of a single junction complex, on average, can be expressed as

(2.14)\begin{align} {\boldsymbol{f}_{\!JC}^{BL/SC}} &= C_{f_{JC}} \left( \boldsymbol{v}^{BL}(\boldsymbol{x}_{JC}) - \boldsymbol{v}^{SC}(\boldsymbol{x}_{JC}) \right) \nonumber\\ &= C_{f_{JC}} \boldsymbol{v}^{BL/SC}(\boldsymbol{x}_{JC}) ={-}C_{f_{JC}} \boldsymbol{v}^{SC/BL}(\boldsymbol{x}_{JC}). \end{align}

The equivalent friction coefficient per unit area of the membrane, denoted as $C_{f}$, is given by

(2.15)\begin{equation} \int_{\triangle S} C_{f} \boldsymbol{v}^{BL/SC}(\boldsymbol{x})\,\mathrm{d}S = \sum_{\boldsymbol{x}_{JC} \in \triangle S} C_{f_{JC}} \boldsymbol{v}^{BL/SC}(\boldsymbol{x}_{JC}) . \end{equation}

For a small patch of surface $\triangle S$, where the velocity variations can be ignored (i.e. $\forall \boldsymbol {x} \in \triangle S$, $\boldsymbol {v}^{BL/SC}(\boldsymbol {x}) = \boldsymbol {v}^{BL/SC}(\boldsymbol {x}_{JC}) \approx {\rm const}.$), we can simplify the expression to

(2.16)\begin{equation} C_{f} \triangle S = C_{f_{JC}} N_{JC} , \end{equation}

where $N_{JC}$ is the number of junction complexes within the patch. We can also introduce the areal density of the junction complex, denoted as $\rho _{JC}$, and obtain

(2.17)\begin{equation} C_{f} = \frac{N_{JC}}{\triangle S} C_{f_{JC}} = \rho_{JC} C_{f_{JC}}. \end{equation}

For our study, we adopted $C_f = 144\ {\rm pN}\ {\rm s}\ \mathrm {\mu }{\rm m}^{-3}$ (Peng et al. Reference Peng, Asaro and Zhu2011). The sensitivity of the double-layer modelling strategies to this parameter is discussed in § 5.

2.7. Cytoskeleton kinematics

The velocity differential between the cytoskeleton and the lipid bilayer is defined by the expression

(2.18)\begin{equation} \boldsymbol{v}^{SC/BL}(\boldsymbol{x}) = \boldsymbol{v}^{SC}(\boldsymbol{x}) - \boldsymbol{v}^{BL}(\boldsymbol{x}) ={-} \frac{1}{C_f} \boldsymbol{f}^{BL/SC}(\boldsymbol{x}) = \frac{1}{C_f} \boldsymbol{\mathsf{I}}^S \boldsymbol{f}^{SC}(\boldsymbol{x}). \end{equation}

To obtain the weak formulation of the problem, we employ the weighted residual method with $\delta \boldsymbol {x}$ as the test function, leading to the expression

(2.19)\begin{equation} \int_S \boldsymbol{v}^{SC/BL}(\boldsymbol{x}) \boldsymbol{\cdot} \delta \boldsymbol{x}\,\mathrm{d}S = \frac{1}{C_f} \int_S \boldsymbol{\mathsf{I}}^S \boldsymbol{f}^{SC}(\boldsymbol{x}) \boldsymbol{\cdot} \delta \boldsymbol{x}\,\mathrm{d}S. \end{equation}

The first member of this equation can be expressed as a function of the mass matrix in terms of the Loop interpolation functions, as given in Boedec et al. (Reference Boedec, Leonetti and Jaeger2017). Notably, the second term in the equation, which represents the internal force density of the cytoskeleton, is similar to the corresponding expression used to determine the elastic membrane forces of the RBC.

3. Comparison of RBC modelling strategies on extensional flow

Laser tweezer stretching is a widely used method to measure the mechanical properties of cells. However, numerical studies have two major shortcomings. Firstly, there is variability in how the opposing forces are applied in simulations (Sigüenza et al. Reference Sigüenza, Mendez and Nicoud2017). Secondly, it is not representative of the stretch that an RBC may experience in a flow since it does not consider the interaction with the surrounding fluid. An alternative approach that avoids these limitations is extensional flow, which has been used for studying vesicles in previous studies (Kantsler, Segre & Steinberg Reference Kantsler, Segre and Steinberg2007, Reference Kantsler, Segre and Steinberg2008; Zhao & Shaqfeh Reference Zhao and Shaqfeh2013; Narsimhan, Spann & Shaqfeh Reference Narsimhan, Spann and Shaqfeh2014; Dahl et al. Reference Dahl, Narsimhan, Gouveia, Kumar, Shaqfeh and Muller2016). In this method, the flow is defined by a single parameter, $\dot {\epsilon }$, which represents the stretching rate. In a Cartesian coordinate system, the velocity component in the direction of stretching ($z$ coordinate) is given by $v_z = \dot {\epsilon } z$. In its planar version, the other two velocity components are $v_x = - \dot {\epsilon } x$ and $v_y = 0$. In its axisymmetric version, they are written as $v_x = - \dot {\epsilon } x /2$ and $v_y = - \dot {\epsilon } y /2$. We considered both configurations, with the axis of symmetry of the RBC along the $y$ axis and a value of $\dot {\epsilon } = 55\ {\rm s}^{-1}$. The viscosity of the surrounding fluid was set at $\eta _{ext} =25\ {\rm mPa}\ {\rm s}$. As the behaviours were found to be similar, we only present the results for the axisymmetric configuration.

We first checked that the double-layer models behave similarly to the single-layer capsule model when the layers are prevented from sliding. This similarity is expected because all models employ the same set of elastic properties, as outlined in table 1. The difference lies in the distribution of these properties. In the capsule–capsule model, the outer capsule mimics the lipid bilayer and assumes responsibility for bending resistance and area conservation. The coefficient of shear elasticity for the outer capsule is considered negligible compared with that of the inner capsule, which solely represents the cytoskeleton properties. The vesicle–capsule model follows a similar approach. However, as the lipid bilayer is now represented by a vesicle, surface incompressibility is treated more rigorously. Consequently, a slightly reduced elongation compared with the other models is observed, reflecting the stiffening effect caused by the surface incompressibility constraint.

When sliding is allowed, the double-layer models deviate significantly from the single-layer capsule model, with a much higher elongation. In dimensionless time $t^*$, the capsule model reaches a state of permanent deformation at $t^*=3$, while the double-layer models continue to stretch out to $t^*=6$. Despite the slightly lower elongation of the vesicle–capsule model due to its better consideration of surface incompressibility, it remains close to that of the capsule–capsule model. The evolution of the RBC's shape is illustrated by the cross-sectional profiles in the three planes of symmetry in figures 2, 3 and 4 at $t^*=0$ (initial biconcave shape), $1$, $3$ and $6$. In all models, the dimple is reduced under stretching, with complete flattening observed in the capsule–capsule model, consistent with its higher elongation.

Figure 2.

Shape evolution in the ($x^*, y^*$) plane perpendicular to the stretching direction for the axisymmetric extensional flow (capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ), for $t^* = 0$ ((a) initial shape – plotted in hereafter), $1$ (b), $3$ (c) and $6$ (d).

Figure 3.

Shape evolution in the ($x^*,z^*$) stretch plane for the axisymmetric extensional flow (capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ), for $t^* = 0$ ((a) initial shape), $1$ (b), $3$ (c) and $6$ (d).

Figure 4.

Shape evolution in the ($y^*,z^*$) stretch plane for the axisymmetric extensional flow (capsule $=$ , vesicle–capsule $=$ , capsule–capsule = ), for $t^* = 0$ ((a) initial shape), $1$ (b), $3$ (c) and $6$ (d).

The mechanical properties of all models are identical, resulting in similar behaviour when sliding is prohibited in the double-layer strategies. However, when sliding is allowed, the observed variation in behaviour must be explained by this additional degree of freedom. To quantify the intensity of sliding, figure 5 displays the dimensionless velocity difference between the cytoskeleton and lipid bilayer, which shows that sliding is maximal in the initial stages and disappears completely when the stationary state is reached. The sliding velocity in the ($x^*, y^*$) plane as a function of $x^*$ (figure 5a) and in the ($y^*,z^*$) plane as a function of $z^*$ (figure 5b) is given for $t^*=1$, $2$ and $6$. The curves of sliding velocity as a function of $z^*$ in the ($x^*,z^*$) plane (not shown here) are similar to those in the ($y^*,z^*$) plane, the intensity is slightly lower, with a maximum of approximately $0.04$ instead of $0.06$, and the sliding is almost identical for both models at $t^*=1$.

Figure 5.

Decrease of the sliding velocity $\delta {\boldsymbol {v}}^*=\boldsymbol {v}^{*BL} - \boldsymbol {v}^{*SC}$ for $t^*=1$, 2 and $6$ for the axisymmetric extensional flow (vesicle–capsule $=$ , capsule–capsule $=$ ), as a function of $x^*$ in the ($x^*,y^*$) plane (a) and of $z^*$ in the ($y^*,z^*$) plane (b).

One explanation for the observed behaviour concerns the degree to which the surface incompressibility of the lipid bilayer is transmitted to the cytoskeleton. In the capsule model, the transmission is total since it is not possible to uncouple the cytoskeleton from the lipid bilayer. This is reflected by a value of the $C^{SK}$ coefficient of the Skalak model as large as possible, i.e. 80 in our simulations. On the other hand, in the double-layer models, the cytoskeleton can relax this constraint, thanks to the possibility of sliding. For these two models, $C^{SK}=2$ for the capsule which represents the cytoskeleton. Figure 6 compares the local relative area variation for the cytoskeleton between the three models. The local area variation is normalised by the total surface area of the RBC at rest, i.e. at the initial time. Although the variations are stronger and very comparable for the vesicle–capsule and capsule–capsule models, this explanation alone does not account for the observed behaviour.

Figure 6.

Relative cytoskeleton area change (represented by the colour code between $-0.25$ in blue and $0.25$ in red) for the axisymmetric extensional flow at $t^* = 6$. $C^{SK} = 80$ for the capsule model (a) and $C^{SK} = 2$ for the capsule representing the cytoskeleton in the vesicle–capsule model (b) and capsule–capsule model (c).

The stress relaxation in the cytoskeleton offered by the possibility of sliding is more general, as it is not just about the local variation in the area. Instead, the cytoskeleton can better manage the whole deformation state imposed on it by the displacement of the RBC surface to which it is subject. The possibility of sliding gives it complete freedom to optimise its elastic stress state by the tangential movement to the surface. Moreover, there is a flow amplification effect, since the intensity of the velocity component according to the direction of the flow increases linearly in $|z^{*}|$. By reducing its deformation energy, the RBC can stretch further, and its two tips venture into regions where the intensity of the stretching velocity is greater, in accordance with the linear increase of the latter.

Figure 7 provides validation for this scenario, with a comparison of the elastic strain energy evolution for the whole RBC membrane (lipid bilayer $+$ cytoskeleton). The graph shows the deformation energy as a function of elongation, which is characterised by the position of the RBC tip in the stretching direction, denoted as $z^*_{max}$. This position is a nonlinear function of time, and the maximum elongation is reached when $z^*_{max}$ stabilises. When sliding is not allowed, the three evolution curves coincide, as seen in the blue curve obtained for the capsule model. In contrast, when sliding is allowed, the curves differ from each other, with the deformation energy being highest for the capsule model and lowest for the capsule–capsule model. The evolution curve for the vesicle–capsule model falls within the envelope formed by the other two curves, illustrating the correlation between the growth of deformation energy and the intensity of elongation. The final $z^*_{max}$ value for the capsule model is lower than 1.7, while it is higher for the vesicle–capsule and capsule–capsule models, reaching almost 1.9 for the latter. The bending contribution to the deformation energy is also plotted on the graph, but its influence is found to be negligible, with the three dashed lines being nearly identical.

Figure 7.

Elastic deformation energy ($\mathcal {W}^{BL}+\mathcal {W}^{SC}$) as a function of the RBC elongation, given by the normalised position of the RBC's tip $z^*_{max}$ at normalised time $t^*$ (capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ). Solid lines: total elastic energy (lipid bilayer $+$ cytoskeleton). Dashed lines: bending energy contribution ($\mathcal {W}^{H}$).

It is interesting to note that the three curves of the total deformation energy all show a minimum for the same elongation value of $z^*_{max} \approx 1.1$, which is related to the chosen reference shape for the cytoskeleton, namely quasi-spherical. The deformation energy of the cytoskeleton decreases as the RBC comes closer to its reference shape, resulting in a minimum value for the deformation energy. If the discocyte reference shape is chosen instead, the curves no longer show a minimum. However, apart from this point, the evolution of the deformation energy is similar for both reference shapes, and the conclusion drawn from them is identical. The elastic strain energy curves are also found to be very close when a spontaneous curvature is considered for the quasi-spherical reference shape. Overall, these findings provide important insights into the relationship between deformation energy and elongation for different RBC models.

4. Comparison of RBC modelling strategies for shear flow

The simple shear flow is the most commonly used configuration for characterising the dynamics of a flowing RBC. By orienting the reference frame such that the $x$ axis corresponds to the direction of flow and the $y$ axis corresponds to that of the velocity gradient, the velocity field can be written as $\boldsymbol {v} = v(y) \boldsymbol {e}_x$, with $v(y) = \dot {\gamma }y$, as shown in figure 8. The $(x, y)$ plane corresponds to the shear plane, and the $z$ axis corresponds to the vorticity axis. The intensity of the shear rate $\dot {\gamma }$ is the single operating parameter that characterises this plane flow. For a suspended particle, such as an RBC, $\dot {\gamma }$ is a hydrodynamic forcing characterised by a viscous stress $\tau = \eta _{ext}\dot {\gamma }$, where $\eta _{ext}$ is the viscosity of the suspending fluid.

Figure 8.

Schematic representation of the shear flow configuration (($x, y, z$) frame) with an RBC (($X, Y, Z$) frame). The angle $\phi$ is the orbit angle, i.e. the angle between the vorticity axis of the flow ($z$ axis) and the symmetry axis of the particle ($Z$ axis). The angle $\theta$ is the inclination angle between the direction of the velocity gradient ($y$ axis) and the projection in the shear plane of the particle's axis of symmetry ($Z$ axis).

Characterising the dynamics of the RBC requires considering its properties as a soft object, in addition to the characteristics of the flow. The first element determining its rigidity is the viscosity contrast $\lambda = \eta _{int}/\eta _{ext}$. Under physiological conditions, $\lambda$ is typically greater than unity because the viscosity of the cytosol of a healthy RBC can vary from $6$ to 20 mPa s (Mohandas & Gallagher Reference Mohandas and Gallagher2008; Williams & Morris Reference Williams and Morris2009), whereas that of the plasma is only approximately $1.5$ mPa s (between $1$ and $1.3$ at body temperature Késmárky et al. Reference Késmárky, Kenyeres, Rábai and Tóth2008). The value of $6$ mPa s is generally considered to be characteristic of the cytosol viscosity of a young and healthy RBC at body temperature. All simulations were carried out with $\eta _{int} = 10\ {\rm mPa}\ {\rm s}$, the value at room temperature. However, most experimental studies have been carried out with much more viscous outside fluids, resulting in characteristic viscosity contrasts below unity.

The second characteristic of the RBC's stiffness is related to the elastic properties of the RBC membrane, via its shear modulus $\mu _s$. Its contribution must be compared with that of the hydrodynamic strength, which gives rise to the introduction of the capillary number $Ca(\dot {\gamma },\mu _s) = \eta _{ext}\dot {\gamma }R/\mu _s = \tau /\tau _{ref}$, the ratio of the viscous stress $\tau = \eta _{ext}\dot {\gamma }$ and $\tau _{ref} = \mu _s/R$, which characterises the intensity of the elastic response of the cytoskeleton. Note that, with our choice of time reference scale $t_{ref} = \eta _{ext}R/\mu _s$, the capillary number may also be defined as $Ca(\dot {\gamma },\mu _s) = \dot {\gamma }t_{ref}$.

The dynamic regimes of RBC, and more generally of a capsule, in shear flow can be represented in the plane $\lambda, Ca(\dot {\gamma },\mu _s)$. For a vesicle, the reference elastic property is the bending modulus $k_b$. The capillary number used therefore is $Ca(\dot {\gamma },k_b) = \eta _{ext}\dot {\gamma }R^3/k_b$, which is related to the capillary number for capsules by the dimensionless Von Kármán number $\mathcal {K} = k_b/R^2\mu _s = (k_b/R^3)/{\tau _{ref}}$, reflecting the relative importance of curvature elasticity (out of plane) vs shear elasticity (in plane). For RBCs, the Von Kármán number is approximately $5\times 10^{-3}$ based on average values from table 1 of Levant & Steinberg (Reference Levant and Steinberg2016).

We use normalised quantities (indicated by a star) and dimensionless input data in our simulations. All surface density of force quantities is normalised by the reference elastic stress $\tau _{ref}$. Hence, the actual values of $\mu _s$ and $k_b$ are not provided. Instead, the former is specified using the capillary number $Ca(\dot {\gamma },\mu _s)$, while the latter is specified using the Von Kármán number $\mathcal {K}$.

Successive experimental studies have continued to enrich the RBC's phase diagram, and a synthesis of the experimental observations is proposed in Minetti et al. (Reference Minetti, Audemar, Podgorski and Coupier2019). Studies often assume the axis of symmetry of the RBC remains in the shear plane. Under these conditions, the RBC's dynamics can be characterised by the inclination angle $\theta$ (see figure 8) and the Taylor deformation parameter $D= (L - B)/(L+B)$, in which $L$ and $B$ are the major and minor axes of the ellipsoid having the same moment of inertia as the RBC. At low values of $\tau\!$, the RBC exhibits a tumbling dynamics. As $\tau$ increases, the axis of symmetry moves closer to the direction of the velocity and the deformation becomes more intense, resulting in the transition from tumbling to tank-treading dynamics. The inclination angle and the Taylor parameter take constant values in the tank-treading regime, with periodic variations around their mean values (swinging) but with decreasing amplitude at higher $\tau\!$. However, experimental studies indicate that an increase in the shear rate is more likely to cause the axis of symmetry to drift out of the shear plane, leading to a rolling dynamics and the appearance of new dynamic regimes showing RBCs with stomatocyte and then multilobe shapes (Mauer et al. Reference Mauer, Mendez, Lanotte, Nicoud, Abkarian, Gompper and Fedosov2018). The critical transition capillary number curve $Ca_c(\dot {\gamma },\mu _s)$ reaches its minimum when $\eta _{ext}$ is greater than 20 mPa s (Fischer & Korzeniewski Reference Fischer and Korzeniewski2013), making it easier to reach the tank-treading regime at high values of $\eta _{ext}$. As we move towards physiological values of $\eta _{ext}$, i.e. at $\lambda > 1$, the curve $Ca_c(\dot {\gamma },\mu _s)$ sharply increases and the transition can only be reached by imposing very high shear rates.

To simplify our analysis, we have assumed that the axis of symmetry of the RBC remains in the shear plane. Although this assumption may not hold for all transition routes, it is sufficient for our purpose, which is to compare modelling strategies at representative points of the RBC dynamics. Moreover, as demonstrated in Levant & Steinberg (Reference Levant and Steinberg2016), the relationship between simple shear flow and plane linear flow allows for this simplified configuration. Another reason to restrict to this simpler configuration is that the characteristic drift time of the RBC symmetry axis out of the shear plane can be large relative to the time scale of RBC dynamics for large $\eta _{ext}$ and small $\dot {\gamma }$ (Levant & Steinberg Reference Levant and Steinberg2016). In addition, under physiological conditions, stable orbits widen as the shear rate decreases, and all orbits become stable below a certain value of $\tau\!$, which is approximately $10^{-2}$ Pa (as shown in figure 14(c) in Minetti et al. Reference Minetti, Audemar, Podgorski and Coupier2019). This trend is consistent at both low and high $\eta _{ext}$ values, as long as we consider equivalent viscous stress. Therefore, exploring the low $\tau$ region is easier when $\eta _{ext}$ is small, as in physiological conditions.

To comprehensively explore the diverse dynamic regimes of an RBC in shear flow, our comparative study focuses on four key points ($\lambda, Ca(\dot {\gamma },\mu _s)$) of the phase diagram. The first point corresponds to the tumbling regime, which is a stable regime for an RBC with its axis of symmetry in the flow plane under physiological conditions. The second point corresponds to the tank-treading regime, which is observed at $\eta _{ext}$ values greater than 20 mPa s. We consider both possible orientations of the axis of symmetry of the RBC, aligned with the axis of vorticity or in the shear plane. The third point is located in the intermittency region, where both tumbling and tank-treading dynamics coexist. Finally, the last point is selected under physiological conditions but with a high shear rate to compare the modelling strategies in the regime of very high deformations.

4.1. Effect of modelling strategy on tumbling dynamics

Our first investigation point considered a viscous stress of 0.01 Pa and an external viscosity of 1.5 mPa s, corresponding to the characteristic dimensionless numbers $\lambda = 6.67$ and $Ca(\dot {\gamma },\mu _s) = 5 \times 10^{-3}$. We found that the inclination angle of the RBC (figure 9a) undergoes a tumbling motion, with the spinning frequency around the vorticity axis varying based on the model used. Interestingly, the capsule model exhibited the highest frequency, while the vesicle model's rotation frequency was almost half that of the other models.

Figure 9.

(a) Time evolution of the inclination angles $\theta /{\rm \pi}$ in tumbling regime at $\lambda = 6.67, Ca(\dot {\gamma },\mu _s) = 5\times 10^{-3}$. Points: simulations (vesicle $=$ black, capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ). Solid lines: Jeffery's theory (second equation of (4.1a,b)) for fitted $r$ on frequency criteria ($r = 2.25$ for capsule, $r = 2.35$ for vesicle–capsule, $r = 2.46$ for capsule–capsule and $r = 4.15$ for vesicle). (b) Time evolution of the deformation parameter $D$.

Additionally, the Taylor deformation parameter $D$ (figure 9b) highlighted a direct relationship between an object's stiffness and its rotation frequency, with the effect being particularly amplified in the vesicle model, which lacks shear elasticity. Further analysis showed that the stiffening effect of the surface incompressibility constraint was clearly demonstrated in the comparison between the capsule–capsule and vesicle–capsule models.

Since the RBC's deformation is small, this dynamic regime is close to that of a rigid particle, and Jeffery's theory (Jeffery Reference Jeffery1922)

(4.1a,b)\begin{equation} r\tan\phi = \frac{C^{orbit}}{\sqrt{r^{{-}2}\cos^2\theta+\sin^2\theta}} , \quad r\tan\theta = \tan \frac{\dot{\gamma}t}{r+r^{{-}1}}, \end{equation}

holds well, where $C^{orbit} = r\tan \phi _0$ is the orbit parameter and $r$ is the particle's aspect ratio. The prolate and oblate shapes are characterised by $r < 1$ and $r > 1$, respectively, with $r$ appearing in the equation $X^2 + Y^2 + r^2Z^2 = 1$ for the surface of an ellipsoidal object. It is important to note that, in this theory, all orbits are probable and determined solely by the initial value $\phi _0$ of the angle $\phi$. To verify Jeffery's theory, we superimposed the curves obtained for each modelling strategy using the second equation of (4.1a,b) as solid lines (figure 9a), with the aspect ratio $r$ adjusted to reproduce the oscillation frequency of the model.

4.2. Effect of modelling strategy on tank-treading dynamics

In the second point, we investigated the effects of increased external viscosity ($\eta _{ext} = 25\ {\rm mPa}\ {\rm s}$) and high shear rate ($\tau = 1.06\ {\rm Pa}$) on RBCs, characterised by the dimensionless numbers $\lambda = 0.4$ and $Ca(\dot {\gamma },\mu _s) = 0.5$. When the RBC axis of symmetry is initially aligned with the axis of vorticity, the transient phase during which it adapts its shape is longer. Figures 10 and 11 show the shape evolution during this phase in the $(x^*,y^*)$ and $(z^*,y^*)$ planes, respectively. All models follow a similar transition pattern with equivalent times, except for the vesicle model, which undergoes a large deformation and then lengthens. However, a difference in behaviour can be observed between single-layer and double-layer strategies. While the capsule model still shows a curvature inversion at $y^* = 0$ in the $(y^*, z^*)$ plane at $t^* = 11$, the double-layer models do not. Figure 12 provides a 3-D view at that time.

Figure 10.

Shape evolutions ($(x^*,y^*)$ plane sectional drawing) in tank-treading regime at $\lambda = 0.4, Ca(\dot {\gamma },\mu _s) = 0.5$ in the transient phase at $t^* = 0$ (a), $t^* = 6$ (b), $t^* = 11$ (c) and $t^* = 15$ (d), for the case when the RBC's symmetry axis is initially aligned with the axis of vorticity (vesicle $=$ black, capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ).

Figure 11.

Shape evolutions ($(y^*,z^*)$ plane sectional drawing) in tank-treading regime at $\lambda = 0.4, Ca(\dot {\gamma },\mu _s) = 0.5$ in the transient phase at $t^* = 0$ (a), $t^* = 6$ (b), $t^* = 11$ (c) and $t^* = 15$ (d), for the case when the RBC's symmetry axis is initially aligned with the axis of vorticity (vesicle $=$ black, capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ).

Figure 12.

Shapes (3-D view) in tank-treading regime at $\lambda = 0.4, Ca(\dot {\gamma },\mu _s) = 0.5$ in the transient phase at $t^* = 11$, for the case when the RBC's symmetry axis is initially aligned with the axis of vorticity. (a) Vesicle, (b) capsule, (c) vesicle–capsule, (d) capsule–capsule.

When the RBC's axis of symmetry is initially in the shear plane, it is already well oriented relative to its final steady state. The transient phase of large deformations is, correspondingly, reduced. The shape evolution during this phase in the $(x^*,y^*)$ and $(z^*,y^*)$ planes is shown in figures 13 and 14, respectively. Figure 15 at $t^* = 15$ provides a 3-D view.

Figure 13.

Shape evolutions ($(x^*,y^*)$ plane sectional drawing) in tank-treading regime at $\lambda = 0.4, Ca(\dot {\gamma },\mu _s) = 0.5$ in the transition phase at $t^* = 0$ (a), $t^* = 6$ (b), $t^* = 11$ (c) and $t^* = 15$ (d), when the RBC's symmetry axis is in the shear plane at the start (vesicle $=$ black, capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ).

Figure 14.

Shape evolutions ($(y^*,z^*)$ plane sectional drawing) in tank-treading regime at $\lambda = 0.4, Ca(\dot {\gamma },\mu _s) = 0.5$ in the transition phase at $t^* = 0$ (a), $t^* = 6$ (b), $t^* = 11$ (c) and $t^* = 15$ (d), when the RBC's symmetry axis is in the shear plane at the start (vesicle $=$ black, capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ).

Figure 15.

Shapes (3-D view) in tank-treading regime at $\lambda = 0.4, Ca(\dot {\gamma },\mu _s) = 0.5$ in the transition phase at $t^* = 15$, when the RBC's symmetry axis is in the shear plane at the start (a) vesicle, (b) capsule, (c) vesicle–capsule, (d) capsule–capsule.

After the transient phase, the dynamics is identical, regardless of the initial orientation of the symmetry axis. The evolution curves in figure 16 are indistinguishable from those obtained when the axis of symmetry is aligned with the axis of vorticity, with a normalised time shift of $\triangle t^* =11.6$. Thus, generalisation to any initial orientation is highly likely.

Figure 16.

Time evolution of the inclination angle $\theta /{\rm \pi}$ (a) and the deformation parameter $D$ (b) in tank-treading regime at $\lambda = 0.4, Ca(\dot {\gamma },\mu _s) = 0.5$, when the RBC's symmetry axis is initially in the shear plane (vesicle $=$ black, capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ).

Once the steady state is reached, the inclination angle's evolution (figure 16a) reveals that the swinging is most pronounced for the capsule model and least pronounced for the vesicle–capsule model. The capsule–capsule model is in between. As expected, the vesicle model exhibits pure tank-treading motion, characterised by a constant value of the inclination angle.

The evolution of deformation via the Taylor parameter $D$ (figure 16b) remains the most effective way to distinguish between the different models. As expected, the vesicle model has the most intense deformation and is completely separate from the other models. The double-layer models produce mean deformations of the same order, with an almost sinusoidal regularity of the Taylor parameter evolution. However, the stiffness provided by the surface incompressibility is reflected in an oscillation amplitude half as large for the vesicle–capsule model as for the capsule–capsule model. Although the mean deformation intensity is lower for the capsule model, its oscillation amplitude is two to three times greater than that of the capsule–capsule model, with the oscillations appearing to be much less symmetrical. We also note that there is an apparent increase in the rotation frequency when the lipid bilayer is modelled as an incompressible fluid film rather than a solid shell. However, it is difficult to conclude whether this increase is related to the consideration of the fluid nature of the lipid bilayer or to the more rigorous treatment of the surface incompressibility constraint.

Figure 17 compares the shape evolutions in the established regime with the vesicle model excluded due to its much larger deformation. The double-layer models produce very similar shapes that are relatively stable, while the capsule model undergoes strong shape variations. In cross-section in the shear plane, the evolution of the capsule model periodically changes from an elliptical to an S-shape (breathing phenomenon).

Figure 17.

Shape evolution in the established tank-treading regime at $\lambda = 0.4, Ca(\dot {\gamma },\mu _s) = 0.5$, when the RBC's symmetry axis is initially in the shear plane. Left panel: $(x^*,y^*)$ plane sectional drawing at $t^* = 7$ (a), $t^* = 15$ (b), $t^* = 20$ (c) and $t^* = 33$ (d) (capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ). Right panel: corresponding 3-D views for capsule (left), vesicle–capsule (middle) and capsule–capsule (right).

A legitimate question is whether the transitional phase in the case where the axis of symmetry is aligned with the axis of vorticity is only a shape adaptation or whether it is accompanied by a general pivoting of the cytoskeleton. If the answer is negative, the final state reached cannot be considered completely identical to the one reached starting from the other orientation, at least from an energetic point of view. To answer this question, we tracked the $y^*$ coordinate of two markers, one initially located on the dimple in the centre of one of the two faces of the RBC and the other on the periphery. The answer is clearly negative: when the axis of symmetry is initially in the shear plane, the dimple undergoes the tank-treading movement, whereas when the axis of symmetry is initially aligned with the axis of vorticity, the marker at the periphery performs the rotation.

4.3. Effect of modelling strategy on the transition to tank treading

Our third point of investigation lies at the upper limit of the intermittency region. For an external viscosity $\eta _{ext}=24\ {\rm mPa}\ {\rm s}$, Fischer & Korzeniewski (Reference Fischer and Korzeniewski2013) identified a critical shear rate of 10 s$^{-1}$ for the transition, corresponding to a critical viscous stress of 0.24 Pa. The corresponding dimensionless numbers are $\lambda = 0.417$ and $Ca(\dot {\gamma },\mu _s) = 0.113$.

Since this point is located in the transition zone, it is expected that the models’ behaviour is particularly sensitive to the modelling strategy. This sensitivity is clearly demonstrated in figure 18. The vesicle model exhibits pure tank treading with a constant inclination angle. The capsule model, on the other hand, stays in the tumbling regime. The double-layer models initially exhibit tumbling, but after one or two periods, they adopt a permanent tank-treading dynamic. Their inclination angle oscillates around a mean value similar to that of the vesicle model, with a slightly larger amplitude for the capsule–capsule model.

Figure 18.

Time evolution of the inclination angle $\theta /{\rm \pi}$ (a) and the deformation parameter $D$ (b) in the tumbling to tank-treading transition region at $\lambda = 0.417, Ca(\dot {\gamma },\mu _s) = 0.113$ (vesicle $=$ black, capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ).

4.4. Effect of modelling strategy in very high deformation regime

Mauer et al. (Reference Mauer, Mendez, Lanotte, Nicoud, Abkarian, Gompper and Fedosov2018) reported a wide range of RBC shapes and dynamics observed in their microfluidic experiments and through simulations using two distinct techniques. Notably, the multilobe regime represents conditions closely resembling the physiological environment, characterised by an external viscosity of $\eta _{ext} = 1.5\ {\rm mPa}\ {\rm s}$, at very high shear rates with $\tau = 3.19$ Pa. This corresponds to dimensionless numbers $\lambda = 6.67$ and $Ca(\dot {\gamma },\mu _s) = 1.5$. We use these simulation parameters to align with the phase diagram presented by Mauer et al. (Reference Mauer, Mendez, Lanotte, Nicoud, Abkarian, Gompper and Fedosov2018), pinpointing a specific point within this multilobe regime. As anticipated, our simulations effectively yielded multilobe shapes, indicating a strong alignment between our modelling approaches and experimental observations. It is worth noting that multilobe shapes have also been observed in a study using a two-dimensional (2-D) vesicle model (Abbasi et al. Reference Abbasi, Farutin, Nait-Ouhra, Ez-Zahraouy, Benyoussef and Misbah2022), suggesting that the cytoskeleton may not be necessary for the multilobe manifestation.

The evolution of the inclination angle and the Taylor deformation parameter is shown in figure 19. The deformation parameter varies between a minimum $D_{min}$ and a maximum $D_{max}$, resulting in different shapes depending on the modelling strategy.

Figure 19.

Time evolution of the inclination angle $\theta /{\rm \pi}$ (a) and the deformation parameter $D$ (b) in the multilobe region at $\lambda = 6.67, Ca(\dot {\gamma },\mu _s) = 1.5$ (capsule $=$ , vesicle–capsule $=$ , capsule–capsule $=$ ). Inserted in (b) are $x$ and $y$ shape views for the capsule and vesicle–capsule models when $D$ is minimal and maximal.

4.5. Comparison with experimental data

In §§ 4.1–4.4, we intentionally selected four points within the phase diagram of an RBC in shear flow. These points were chosen to comprehensively cover the spectrum of the RBC dynamics. Across each of these chosen points, our simulations consistently replicated the anticipated shapes and dynamic responses, including the characteristic multilobe shape. These consistent agreements between our simulation results and experimental observations serve as compelling validation, affirming the effectiveness and reliability of our modelling approaches.

In this subsection, we directly compare our numerical results with existing experimental data. Specifically, we analyse the relationship between the tank-treading frequency ($\,f$) and the shear rate ($\dot {\gamma }$). Fischer (Reference Fischer2007) previously reported that the frequency scales with shear rate as $f\sim \dot {\gamma }^{\beta }$, with the scaling exponent $\beta$ in the range of 0.85 to 0.95.

Employing the same viscosity value for the suspending medium ($\eta _{ext} = 28.9\ {\rm mPa}\ {\rm s}$) as reported in Fischer (Reference Fischer2007), which leads to $\lambda = 0.346$, we conducted a series of simulations using both single-layer and double-layer models to investigate the tank-treading frequency as a function of shear rate. Figure 20 displays a comparison between our simulation results and experimental data. The results show that both single-layer and double-layer models indeed exhibit a power-law relationship with $\beta$ around 0.92. We note that our numerical results are also consistent with the results obtained by Peng et al. (Reference Peng, Li, Pivkin, Dao, Karniadakis and Suresh2013), where they reported a value of approximately 0.91 for both their one-component and two-component models. In our simulations, we assumed a quasi-spherical reference shape with non-zero spontaneous curvature. We also explored scenarios with zero spontaneous curvature, resulting in negligible deviations from the primary results.

Figure 20.

The frequency ($\,f$) of the tank-treading motion of an RBC in shear flow plotted against the shear rate ($\dot {\gamma }$). Simulated results using various models are denoted by filled symbols, while experimental data from Fischer (Reference Fischer2007) are represented by open circles. An exponential fit to the numerical results ($\,f\sim \dot {\gamma }^{0.92}$) is shown as a dashed line, and a linear fit ($\,f\sim \dot {\gamma }$) is depicted by the solid line.