Finite-Element Method for the Simulation of Lipid Vesicle/Fluid Interactions in a Quasi–Newtonian Fluid Flow

Abstract

We present a computational framework for modeling an inextensible single vesicle driven by the Helfrich force in an incompressible, non-Newtonian extracellular Carreau fluid. The vesicle membrane is captured with a level set strategy. The local inextensibility constraint is relaxed by introducing a penalty which allows computational savings and facilitates implementation. A high-order Galerkin finite element approximation allows accurate calculations of the membrane force with high-order derivatives. The time discretization is based on the double composition of the one-step backward Euler scheme, while the time step size is flexibly controlled using a time integration error estimation. Numerical examples are presented with particular attention paid to the validation and assessment of the model’s relevance in terms of physiological significance. Optimal convergence rates of the time discretization are obtained.

1. Introduction

Motivated by their biomedical applications and their great potential for improving existing therapies, the development of accurate numerical tools to better understand the functioning processes of biomemetic lipid vesicles remains a dynamic and stimulating area of research [1]. They represent an attractif model system for a quantitative and predictive understanding of red blood cells (RBCs), whose dynamics continue to pose a formidable challenge to computational modeling [2]. Blood in the microvasculature is a complex fluid whose physiological behavior is critically dependent on its major cellular components, the RBCs. Although blood flow is generally assumed to be Newtonian fluid in the large vascular system, it exhibits a non-Newtonian shear-thinning rheology that cannot be overlooked at the microscopic scale, due to the complex mixture in the plasma [3,4,5,6,7]. Among several existing constitutive laws, we consider the Carreau’s constitutive law, commonly used for blood flow. That is, the viscosity exponentially decreases with increasing shear rates. A major challenge in complex flows lies in deciphering the fluid/membrane space-time interaction at the individual level [8]. Therefore, the simulation of such biomembranes needs to gain in accuracy. This work is oriented towards the direct numerical simulation of a deformable vesicle in a surrounding non Newtonian flow.

A variety of numerical methods have been developed for the vesicle problem by tracking explicitly (e.g., boundary integral method [9,10,11], immersed boundary method [12], penalty-immersed boundary method [13] lattice Boltzmann method [14], parametric finite element method [15,16], etc.) or implicitly (e.g., level set method [17,18,19,20,21], phase field approach [22,23], isogeometric phase field method [24], combined level set and phase field methods [25], etc.) the vesicle’s membrane. However, most numerical studies have been limited to Newtonian fluid models in the zero Reynolds limit. Nevertheless, an astonishing change in the RBC behaviors under finite Reynolds regimes was first pinpointed in [18,26]. To the best of our knowledge, the most relevant numerical study of the dynamics of a single Newtonian vesicle in a non-Newtonian fluid is presented in  [12], where an Oldroyd-B viscoelastic model was used to model the surrounding fluid of the vesicle.

To model the mechanical properties of the membrane of vesicles (and of biological membranes in general), a physiologically realistic model consists of the minimization of the scalar energy functional known as the Helfrich energy. In this model, the main mode of deformation is bending and the cost in bending is determined by the main curvatures of the membrane. That results in a highly nonlinear membrane force with respect to the shape. The model has been introduced in the early 1970s in separate works by Canham [27], Helfrich [28,29] and Evans [30]. For a phospholipidic membrane $\mathsf{\Gamma }$, the Helfrich energy is expessed with respect to the mean curvature H as follows:

$$E(\mathsf{\Gamma })=\frac{k_b}{2}{\int }_{\mathsf{\Gamma }}{\left(H(\mathsf{\Gamma })\right)}^2\mathrm{d}\mathit{s},$$

(1)

with $k_b\approx {10}^{-20}\mathrm{kg}{\mathrm{m}}^2{\mathrm{s}}^{-2}$ the bending rigidity modulus [31,32]. Moreover, a phospholipid bilayer is purely fluidic interface characterized by the local inextensibility that constraints the surface divergence of the velocity field to vanish on the membrane. Membrane incompressibility plays an important role in the dynamics of flowing vesicles and leads to the preservation of the membrane surface (or perimeter in two-dimensional studies) thanks to Reynolds’ lemma, see [18]. This constraint can be enforced, for example, via an exact Lagrange multiplier localised on the membrane which can be thought of as a position-dependent membrane tension [33] or an elastic force expressed in terms of membrane stretching [17]. The fluid/membrane coupling remains challenging due to the highly nonlinear membrane force with respect to the the shape of the vesicle. Instead of using a body-fitted grid for the membrane description, we employ a level set Eulerian framework for a straightforward representation of the sharp interface with a signed distance function [34,35,36]. Generally speaking, the level set method makes it possible to naturally handle the topological changes of interfaces such as interface merging or breaking. On the contrary, numerical problems related to the mesh distortion can raise when modeling the large displacement of an interface or topological changes using Lagrangian approaches. Special treatments and remeshing procedures are required in such situations, which generally remain difficult. The mechanics of the membrane is introduced through a localized body-force term [19]. The deflation of the vesicles is quantified by a dimensionless geometric parameter: the reduced area

$$v=\frac{(\mathcal{A}/\pi )}{{(\mathcal{P}/2\pi )}^2}⩽1,$$

with $\mathcal{P}$ the perimeter of the vesicle and $\mathcal{A}$ the area of the enclosed subdomain.

We have arranged the remainder of the article as follows. Section 2 presents the mathematical framework for an individual vesicle immersed in a non-Newtonian Carreau fluid and the nonlinear coupled problem. It also introduces the penalty method and details the variational formulation of the problem. Section 3 is devoted to the time discretization scheme and provides the global algorithm used in the numerical simulations. In Section 4, we report sample numerical simulations in the Newtonian and non-Newtonian cases with some validation results to illustrate the accuracy of the method. The conclusion and future research directions are presented in Section 5.

2. Problem Statement

The spatio-temporal deformations of the vesicle are driven by the bending force, actions of fluid forces and boundary conditions, requiring the balance equations of mass and momentum. The fluid-membrane coupling is described through the balance of hydrodynamic stress by the bending response of the membrane. The level set description of the membrane obeys a Hamilton–Jacobi equation.

2.1. Level Set Formulation

Let $T>0$ and for any time $t\in (0,T)$, consider a 2D vesicle $\mathsf{\Gamma }(t)=\partial \mathsf{\Omega }(t)$ immersed in a non Newtonian fluid occupying both extracellular $\mathsf{\Lambda }\backslash \mathsf{\Omega }(t)$ and enclosed $\mathsf{\Omega }(t)$ spaces. Here, $\mathsf{\Lambda }\in {\mathbb{R}}^2$ is a bounded convex polygonal domain with $\mathsf{\Gamma }(t)\cap \partial \mathsf{\Lambda }=\varnothing$. See Figure 1. From now, the explicit dependence of $\mathsf{\Gamma }(t)$ and $\mathsf{\Omega }(t)$ from time will be understood but will be omitted.

The membrane is implicitely described as a zero-level set of a Lipschitz continuous scalar function as follows

$$\mathsf{\Gamma }=\left\{(t,\mathit{x})\in (0,T)\times \mathsf{\Lambda }:\phi (t,\mathit{x})=0\right\}.$$

We assume that the region with negative $\phi (t,.)$ represents the internal domain of the vesicle $\mathsf{\Omega }$. Generally speaking, the main utility of an implicit representation of an interface is that large deformations or changes in topology, such as merging or pinching, which frequently occur when dealing, for example, with bubbles and drops, can be treated in a natural and direct way without the need for meshing or local refinement of the mesh in the vicinity of the complex interface; Authomatic meshing remains challenging especially in three-dimensional setups. The level set function is initialized as a signed distance function ${\phi }_0$ with respect to the initial position of the vesicle $\mathsf{\Gamma }$. It is given by

$${\phi }_0(\mathit{x})=\left\{\begin{array}{c}{\displaystyle \underset{{\displaystyle \mathit{y}\in \mathsf{\Gamma }(0)}}{inf}\parallel \mathit{y}-\mathit{x}\parallel ,\mathrm{if}\mathit{x}\notin \mathsf{\Omega }(0)}\\ {\displaystyle -\underset{{\displaystyle \mathit{y}\in \mathsf{\Gamma }(0)}}{inf}\parallel \mathit{y}-\mathit{x}\parallel ,\mathrm{otherwise}.}\end{array}\right.$$

Let $\mathbf{I}$ be the identity tensor. Geometric quantities, such that the unit normal outward vector $\mathbf{n}=\nabla \phi /|\nabla \phi |$, the mean curvature $H=\mathrm{div}\mathbf{n}$, the surface projector ${\mathit{\pi }}_{\mathsf{\Gamma }}=\mathrm{I}-\mathbf{n}\otimes \mathbf{n}$, as well as the surface gradient ${\mathbf{\nabla }}_s={\mathit{\pi }}_{\mathsf{\Gamma }}\mathbf{\nabla }$, the surface divergence ${\nabla }_s·=tr({\mathbf{\nabla }}_s)$ and the surface Laplacian ${\Delta }_s$ operators are encoded in terms of $\phi$. These quantities are consequently extended to the entire computational domain $\mathsf{\Lambda }$. The Helfrich force requires a fourth-order derivative of $\phi$, yielding a numerically stiff problem with severe time-step limitations needed for stability [33,37].

2.2. Fluid Constitutive Equations

Let $\mathit{u}$ and p be the fluid velocity and pressure fields. Consider a quasi Newtonian constitutive law, where the Cauchy stress

$$\mathit{\sigma }=2\eta \left({\left|\mathbf{D}(\mathbf{u})\right|}^2\right)\mathbf{D}(\mathbf{u})-p\mathbf{I}$$

is expressed with respect to the rate-of-strain

$$\mathbf{D}(\mathbf{u})=\frac{1}{2}(\mathbf{\nabla }\mathit{u}+\mathbf{\nabla }{\mathit{u}}^T).$$

Let us consider the Euclidean norm of tensors $\left|\mathit{\tau }\right|={(\frac{1}{2}\mathit{\tau }:\mathit{\tau })}^{1/2}$. The shear rate is then given by $\left|2\mathbf{D}(\mathbf{u})\right|$. The shear-rate dependent viscosity follows the Carreau’s law

$$\eta \left(\gamma \right)={\eta }_{\infty }+({\eta }_0-{\eta }_{\infty }){\left(1+\kappa \gamma \right)}^{{\displaystyle \frac{\varsigma -1}{2}}},$$

with $\gamma \in {\mathbb{R}}^{+},{\eta }_0⩾{\eta }_{\infty }>0$; See ([38] Section 2.1, Definition 2.3). Note that we followed the notations in ([38] Definition 2.5), where the viscosity function is expressed with respect to the square of the shear rate $\gamma ={\left|2\mathbf{D}(\mathbf{u})\right|}^2$, instead of the shear rate directly. Blood flow features a shear thinning behavior, that is $\varsigma <1$, which implies that the viscosity decreases with the shear rate. We consider the physiological parameters of human blood (assuming a hematocrit concentration $45\%$ at an ambient temperature ${37}^{\circ }$), obtained by fitting data from normal blood samples in [4]. We assume piecewise constant zero-shear viscosity ${\eta }_0={\mu }_i$ and ${\eta }_0={\mu }_o=0.16\mathrm{Pa}\mathrm{s}$ in the enclosed $\mathsf{\Omega }$ and surrounding $\mathsf{\Lambda }\backslash \mathsf{\Omega }$ subdomains, respectively.  The latter correspond to the viscosities in the Newtonian case $\varsigma =1$. Without loss of generality, we assume equal infinite shear viscosities ${\eta }_{\infty }=0.0035\mathrm{Pa}\mathrm{s}$ of the inner and outer fluids.

A linear shear flow is characterized by the imposition of a no-slip condition with constant opposing velocity $\mathit{u}={\mathit{u}}_b$ on the horizontal boundaries ${\Sigma }_D$ of a square domain, while stress-free boundary conditions are prescribed elsewhere on ${\Sigma }_N=\partial \mathsf{\Lambda }\backslash {\Sigma }_D$ [39].

Although $\phi$ is initialized as a signed distance function with respect to the membrane of the vesicle, the advection of the level set function degenerates its initial signed distance property and the norm of the level set gradient may locally vanish or blow up, resulting in numerical singularities. This issue is worked around by solving an auxiliary redistancing equation which helps recover the original signed distance behavior [34,40]. Only a few redistancing iterations are usually done on a regular basis, as the level set function must be kept close to a signed distance map and not necessarily an exact distance, see, e.g., [37].

Three dimensionless physical parameters arise: The Reynolds number Re compares the strength of inertial forces to the viscosity effect. The bending number Bn relates the force of the imposed flow to the characteristic membrane bending strength. The viscosity contrast $\overline{\mu }={\mu }_i/{\mu }_o$ varies by changing the viscosity of the inner fluid and is relevant for the vesicle’s dynamics. Equal densities are assumed in inner and outer domains as this has no effect on the dynamics of the BRC [41]. We introduce the smoothed Heaviside ${\mathcal{H}}_{\epsilon }(·)$ and Dirac ${\delta }_{\epsilon }(·)$ functions within a banded strip of width $2\epsilon >0$ around $\mathsf{\Gamma }$ for easy approximations of integrals. They regularized functions are given by

$$\begin{array}{ccc}\hfill {\mathcal{H}}_{\epsilon }(\phi )& =& \left\{\begin{array}{cc}0,\hfill & \mathrm{if} \phi <-\epsilon \hfill \\ {\displaystyle \frac{1}{2}\left({\displaystyle 1+\frac{\phi }{\epsilon }+\frac{sin\left({\displaystyle \frac{\pi \phi }{\epsilon }}\right)}{\pi }}\right),}\hfill & \mathrm{if} \left|\phi \right|⩽\epsilon \hfill \\ 1,\hfill & \mathrm{otherwise}\hfill \end{array}\right.,\hfill \\ \hfill {\delta }_{\epsilon }(\phi )& =& {\displaystyle \frac{\mathrm{d}{\mathcal{H}}_{\epsilon }}{\mathrm{d}\phi }(\phi )=\left\{\begin{array}{cc}{\displaystyle \frac{1}{2\epsilon }\left(1+cos\left(\frac{\pi \phi }{\epsilon }\right)\right),}\hfill & \mathrm{if} \left|\phi \right|⩽\epsilon \hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right..}\hfill \end{array}$$

The standard finite element approximation would result in suboptimal convergence due to the level set regularization. See for example [42] for more details on this issue.

From now, all quantities are non-dimensionalized. The same notations will be used for ease of exposition. The normalized regularized viscosity is expressed as

$${\eta }_{\epsilon }(\gamma ,\phi )={\eta }_{\infty }+\left(\left({\mathcal{H}}_{\epsilon }(\phi )+\overline{\mu }\left(1-{\mathcal{H}}_{\epsilon }(\phi )\right)\right)-{\eta }_{\infty }\right){\left(1+\kappa \gamma \right)}^{{\displaystyle \frac{\varsigma -1}{2}}},\mathrm{with}\gamma >0.$$

2.3. Statement of the Fluid-Vesicle Interaction Problem

Assume fluid incompressibility (mass conservation (2c)). Completed with the conditions of continuity and mechanical equilibrium at the interface, the normalized coupled fluid/membrane problem takes the form: Find $\phi$, $\mathit{u}$, p and $\lambda$ such that

$$\begin{array}{ccc}\hfill \frac{\partial \phi }{\partial t}+\mathit{u}.\mathbf{\nabla }\phi & =& 0,(0,T)\times \mathsf{\Lambda },\hfill \end{array}$$

(2a)

$$\begin{array}{ccc}\hfill \mathrm{Re}\left(\frac{\partial \mathit{u}}{\partial t}+\mathit{u}.\mathbf{\nabla }\mathit{u}\right)-\mathbf{\nabla }·\left({\eta }_{\epsilon }\left({\left|2\mathbf{D}(\mathbf{u})\right|}^2,\phi \right)\right)+\mathbf{\nabla }p& =& \mathbf{0},(0,T)\times \mathsf{\Lambda },\hfill \end{array}$$

(2b)

$$\begin{array}{ccc}\hfill \nabla .\mathit{u}& =& 0,(0,T)\times \mathsf{\Lambda },\hfill \end{array}$$

(2c)

$$\begin{array}{ccc}\hfill {\nabla }_s.\mathit{u}& =& 0,(0,T)\times \mathsf{\Gamma }.\hfill \end{array}$$

(2d)

Let $[[·]]$ be the discontinuity jump of a given field across the vesicle membrane $\mathsf{\Gamma }$ in its normal direction $\mathbf{n}$. The system (2a–2d) is endowed with the conditions of continuity and mechanical equilibrium at the interface

$$\begin{array}{ccc}\hfill [[\mathit{u}]]& =& \mathbf{0},\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill [[\mathit{\sigma }\mathbf{n}]]& =& {\mathbf{\nabla }}_s\lambda -\lambda H\mathbf{n}+\frac{1}{\mathrm{Bn}}\left(2{\Delta }_{s}H+\frac{H^3}{2}\right)\mathbf{n}.\hfill \end{array}$$

(4)

The latter is transferred as a source term in the momentum equation. We refer to [18] for more details on the total stress jump across the vesicle and the derivation of the membrane bending force. Moreover, we consider the boundary conditions $\phi ={\phi }_b$ and $\mathit{u}={\mathit{u}}_b$ on $(0,T)\times {\Sigma }_{-}$, and $\mathit{\sigma }·\mathbf{n}=\mathbf{0}$ on $(0,T)\times {\Sigma }_N$. Here ${\Sigma }_{-}(t)=\{\mathit{x}\in \partial \mathsf{\Lambda }:\mathit{u}(t,\mathit{x})·\mathit{\nu }(\mathit{x})<0\}$ is the inflow and $\mathit{\nu }$ is the outward-pointing unit normal to $\mathsf{\Lambda }$, see Figure 1.

2.4. Panalty Approach

For the shear flow setups, we introduce $\mathbb{V}({\mathit{u}}_b)=\{\mathit{v}\in {(H^1(\mathsf{\Lambda }))}^2:\mathit{v}={\mathit{u}}_b \mathrm{on} {\Sigma }_D\}$. Hence, the admissible velocities belong to the constrained space

$$\mathbb{K}({\mathbf{u}}_b)=\left\{\mathbf{v}\in \mathbb{V}({\mathit{u}}_b):\nabla .\mathit{v}=0 \mathrm{in} \mathsf{\Lambda },{\nabla }_s.\mathit{v}=0 \mathrm{in} \mathsf{\Gamma }\right\}.$$

Analogously to the pressure p introduced to impose the incompressibility of the fluid, the zero surface divergence constraint (2d) is imposed by an exact local Lagrange multiplier $\lambda$ on $\mathsf{\Gamma }$ [33]. Targeting reduced computational load and simple implementation using standard fluid solvers, a penalty method relaxes the inextensibility constraint. For ease of presentation, the penalized problem is presented for vanishing inertia. The unsteady case can be written after time discretization, see a detailed presentation for a different problem in [43]. Let us introduce the dimensionless sharp function:

$${\eta }_0^{*}(t,\mathit{x})=\left\{\begin{array}{cc}\hfill \overline{\mu }& \mathrm{when} \mathit{x}\in \mathsf{\Omega }(t),\hfill \\ \hfill 1& \mathrm{otherwise} \mathrm{for} \mathsf{\Lambda }\backslash \mathsf{\Omega }(t).\hfill \end{array}\right.$$

The problem (2) is first witten as a minimization under constraint:

$$(\mathbf{u},\mathsf{\Omega })={\displaystyle \underset{\begin{array}{c}\mathbf{v}\in \mathbb{K}(t,{\mathbf{u}}_b),\\ \omega \subset \mathsf{\Lambda }\end{array}}{\arg \inf }}J(\mathbf{v},\omega ),$$

(5)

with a total energy functional

$$J(\mathbf{v},\omega )={\int }_{\mathsf{\Lambda }}\mathcal{J}\left({\left|2\mathbf{D}(\mathbf{u})\right|}^2\right)\mathrm{d}x+\frac{1}{2\mathrm{Bn}}{\int }_{\partial \omega }{H}^2(\omega )\mathrm{d}s.$$

(6)

Here, the scalar differentiable energy of dissipation reads

$$\mathcal{J}(\gamma )=\frac{{\eta }_{\infty }}{2}\gamma +\frac{({\eta }_0^{*}-{\eta }_{\infty })}{\kappa (\varsigma +1)}\left({(1+\kappa \gamma )}^{{\displaystyle \frac{\varsigma -1}{2}}}-1\right).$$

Note that the Helfrich energy term depends indirectly on the velocity field $\mathbf{u}$, since the function $\phi$ describing the shape, is advected by $\mathbf{u}$. The optimality system results from a saddle-point formulation. To alleviate the inextensibility, we introduce the space of unconstrained incompressibility

$$\tilde{\mathbb{K}}({\mathbf{u}}_b)=\{\mathbf{v}={\mathit{u}}_b:\nabla ·\mathit{v}=0 \mathrm{in} \mathsf{\Lambda }\}$$

For a small penalty parameter $ϵ$, the minimization problem (5) is approached by another minimization problem penalizing (2d):

$$(\mathbf{u},\mathsf{\Omega })=\underset{\begin{array}{c}\mathbf{v}\in \tilde{\mathbb{K}}(t,{\mathbf{u}}_b)\\ \omega \subset \mathsf{\Lambda }\end{array}}{\arg \inf }{J}_{ϵ}(\mathbf{v})=J(\mathit{v})+\frac{1}{ϵ}{\int }_{\partial \omega }{\left({\nabla }_s·\mathit{v}\right)}^2\mathrm{d}s.$$

Using a saddle point formulation, the optimality system produces a penalty term in the momentum equation. The level set advection equation is solved by the Streamline Upwind Petrov-Galerkin method, with a stabilisation term $S(\phi ,\psi )$ introducing diffusion in the streamline direction [19].

Spatial discretization is carried out by finite elements. To reduce spurious oscillations and errors due to multi-space projections for lower degree polynomials and accurately calculate the bending term, higher order finite elements are considered using Lagrange polynomials of degrees $k+1$, k and $k+1$ for $\mathit{u}$, p and $\phi$, respectively, with $k=\{1,2,3\}$. Considering the total derivative term, the variational formulation of the panalized problem writes:

Find $\mathit{u}\in {\mathcal{C}}^0((0,T),\mathbb{V}({\mathit{u}}_b))$, $p\in L^2((0,T),L_0^2(\mathsf{\Lambda }))$ and $\phi \in {\mathcal{C}}^0((0,T),W^{1,\infty }\left(\mathsf{\Lambda }\right)\cap H^1(\mathsf{\Lambda }))$ such that

$$\begin{array}{ccc}& & \mathrm{Re}{\int }_{\mathsf{\Lambda }}\left(\frac{\partial \mathit{u}}{\partial t}+\mathit{u}·\nabla \mathit{u}\right)·\mathit{v}+{\int }_{\mathsf{\Lambda }}2{\eta }_{\epsilon }\left({\left|2\mathbf{D}(\mathbf{u})\right|}^2,\phi \right)\mathbf{D}(\mathbf{u}):\mathbf{D}(\mathbf{v})-{\int }_{\mathsf{\Lambda }}p\nabla .\mathit{v}\hfill \end{array}$$

$$\begin{array}{ccc}& & +\frac{1}{ϵ}{\int }_{\mathsf{\Lambda }}|\nabla \phi |{\delta }_{\epsilon }(\phi )({\nabla }_s.\mathit{u})({\nabla }_s.\mathit{v})-\frac{1}{\mathrm{Bn}}{\int }_{\mathsf{\Lambda }}|\nabla \phi |{\delta }_{\epsilon }(\phi )\left(\mathsf{\Psi }+\frac{H^3}{2}\right)\mathbf{n}·\mathit{v}=0,\hfill \end{array}$$

(7a)

$$\begin{array}{ccc}& & {\int }_{\mathsf{\Lambda }}q\nabla .\mathit{u}=0,\hfill \end{array}$$

(7b)

$$\begin{array}{ccc}& & {\int }_{\mathsf{\Lambda }}|\mathbf{\nabla }\phi |\mathsf{\Psi }\xi +{\int }_{\mathsf{\Lambda }}|\mathbf{\nabla }\phi |\left({\nabla }_{s}H\right)·\left({\nabla }_s\xi \right)=0,\hfill \end{array}$$

(7c)

$$\begin{array}{ccc}& & {\int }_{\mathsf{\Lambda }}\frac{\partial \phi }{\partial t}\psi +{\int }_{\mathsf{\Lambda }}\left(\mathit{u}·\mathbf{\nabla }\phi \right)\psi +{\int }_{\mathsf{\Lambda }}S\left(\phi ,\psi \right)=0,\hfill \end{array}$$

(7d)

$\forall \mathit{v}\in \mathbb{V}(\mathbf{0})$, $\forall q\in L_0^2(\mathsf{\Lambda })$, $\forall \xi \in H^1(\mathsf{\Lambda })$   and   $\forall \psi \in H^1(\mathsf{\Lambda })$.

3. Time Discretization

We consider a first order backward Euler scheme, referred to as BE, usually used for the vesicle’s problem, and rely on the flow composition technique to raise the order of the time integration scheme [44]. An adaptive time-stepping will be introduced afterwards.

3.1. Error Estimation

Henceforth, we focus on the time discretization of the system of first-order differential equations obtained by finite element discretization. It is simply denoted by

$${\mathit{y}}^{\prime }=\mathit{f}(t,\mathit{y}(t)),\mathrm{with}\mathit{y}(0)={\mathit{y}}_0\mathrm{and}t\in \mathbb{R}.$$

(7e)

Assume f is an analytic function. Let $t_n,0⩽n⩽N$ be a partition of the time interval $[0,T]$ with time step $\delta t=t_{n+1}-t_n$. and ${\mathit{y}}_n$ be the approximations of the solution $\mathit{y}(t_n)$. Let ${𝓨}_t:\mathbb{R}\to \mathbb{R}$ be the exact t-flow of the above differential equation so that with ${𝓨}_t({\mathit{y}}_0)=\mathit{y}(t)$. A basic integrator, such that ${\mathit{y}}_{n+1}={\mathit{Y}}_{\delta t}({\mathit{y}}_n)$, is said of order p if the following equation holds

$${𝓨}_{\delta t}({\mathit{y}}_0)-{\mathit{Y}}_{\delta t}({\mathit{y}}_0)=\mathcal{O}\left({\delta t}^{p+1}\right).$$

The order p of a one-step method can be raised by composition ${\mathit{Y}}_{a_1\delta t}\circ \dots \circ {\mathit{Y}}_{a_s\delta t}$ with step sizes $a_1\delta t,\cdots ,a_s\delta t$ to at least an order $p+1$ if the conditions

$$a_1+\cdots +a_s=1\mathrm{and}{a}_1^{p+1}+\cdots +a_s^{p+1}=0$$

are satisfied [44]. No real solutions hold for odd p; This can be circumvented by considering, for instance, an analytic continuation off the real axis by using complex coefficients $a_i$’s [45]. For a thorough description of the composition of multi-step methods and the choice of parameters $a_1,\cdots ,a_s$ for any general order p and step number s, we refer to the book ([44] Section II.4).

In what follows, we limit our approach to the double composition (i.e., we set $s=2$) of the numerical flow ${\mathit{Y}}_{\delta t}$ associated with the backward Euler scheme, that is

$${\mathit{Y}}_{\delta t}({\mathit{y}}_n)={\mathit{y}}_n+\delta t\mathit{f}(t_n,{\mathit{y}}_n).$$

The aforementioned algebraic conditions provide the coefficients $2a_1=2{\overline{a}}_2=1+i$.

Theorem 1.

For $n⩾0$, consider the first order differential equations in $[t_n,t_{n+1}]$. Let us denote ${\mathit{y}}_{n+1,0}={\mathit{y}}_n$, ${\mathit{y}}_{n+1,1}={\mathit{Y}}_{a_1\delta t}\left({\mathit{y}}_{n+1,0}\right)$ and ${\mathit{y}}_{n+1,2}={\mathit{Y}}_{a_2\delta t}\left({\mathit{y}}_{n+1,1}\right)={\mathit{Y}}_{a_2\delta t}\circ {\mathit{Y}}_{a_1\delta t}\left({\mathit{y}}_n\right)$. The composed method defined by the real part ${\mathit{y}}_{n+1}=Re\left({\mathit{y}}_{n+1,2}\right)$ is a numerical method of order 2. Furthermore, the imaginary part, referred to as $Im\left({\mathit{y}}_{n+1,2}\right)$, represents an error estimate of the approximation to within $\mathcal{O}\left({\delta t}^3\right)$.

Proof. 

For $n⩾1$ and $k=1,2$, we shall expand $\mathit{f}(t_n+a_k\delta t,{\mathit{y}}_{n+1,k})$ into a Taylor series around $(t_n,{\mathit{y}}_{n+1,k-1})$. We have

$${\mathit{y}}_{n+1,1}={\mathit{y}}_{n+1,0}+a_1\delta t\mathit{f}\left(t_n,{\mathit{y}}_{n+1,0}\right)$$

and

$${\mathit{y}}_{n+1,2}={\mathit{y}}_{n+1,1}+a_2\delta t\mathit{f}\left(t_n+a_1\delta t,{\mathit{y}}_{n+1,1}\right).$$

Considering ${\mathit{y}}_{n+1,2}$, we expand the right-hand expression into a Taylor series and truncate

$$\begin{array}{cc}\hfill {\mathit{y}}_{n+1,2}& ={\mathit{y}}_{n+1,1}+a_2\delta t\mathit{f}\left(t_n+a_1\delta t,{\mathit{y}}_{n+1,1}\right)\hfill \\ \hfill & ={\mathit{y}}_n+a_1\delta t\mathit{f}(t_n,{\mathit{y}}_n)+a_2\delta t\mathit{f}\left(t_n+a_1\delta t,{\mathit{y}}_n+a_1\delta t\mathit{f}(t_n,{\mathit{y}}_n)\right)\hfill \\ \hfill & ={\mathit{y}}_n+a_1\delta t\mathit{f}(t_n,{\mathit{y}}_n)\hfill \\ \hfill & +a_2\delta t[\mathit{f}(t_n,{\mathit{y}}_n)+a_1\delta t\left({\mathit{f}}_t+\mathit{f}{\mathit{f}}_{\mathit{y}}\right)(t_n,{\mathit{y}}_n).\hfill \\ \hfill & +\frac{1}{2}{a}_1^2\delta t^2\left({\mathit{f}}_{t,t}+2{\mathit{f}}_{t,\mathit{y}}\mathit{f}+{\mathit{f}}_{\mathit{y},\mathit{y}}{\mathit{f}}^2\right)(t_n,{\mathit{y}}_n)+\mathcal{O}\left(\delta t^3\right)]\hfill \\ \hfill & ={\mathit{y}}_n+(a_1+a_2)\delta t\mathit{f}(t_n,{\mathit{y}}_n)+(a_1{a}_2)\delta t^2\left({\mathit{f}}_t+{\mathit{f}}_{\mathit{y}}\mathit{f}\right)(t_n,{\mathit{y}}_n)\hfill \\ \hfill & +\frac{1}{2}{a}_2{a}_1^2\delta t^2\left({\mathit{f}}_{t,t}+2{\mathit{f}}_{t,\mathit{y}}\mathit{f}+{\mathit{f}}_{\mathit{y},\mathit{y}}{\mathit{f}}^2\right)(t_n,{\mathit{y}}_n)+a_2\delta t\mathcal{O}\left(\delta t^3\right).\hfill \end{array}$$

The subscript refers to the partial derivatives with respect to the given variable. Using the expression of ${\mathit{y}}^{(2)}$, we get

$$\begin{array}{cc}\hfill {\mathit{y}}_{n+1,2}& ={\mathit{y}}_n+(a_1+a_2)\delta t{\mathit{y}}^{\prime }(t_n)+(a_1{a}_2)\delta t^2{\mathit{y}}^{″}(t_n)\hfill \\ \hfill & +\frac{1}{2}{a}_2{a}_1^2\delta t^2\left({\mathit{f}}_{t,t}+2{\mathit{f}}_{t,\mathit{y}}\mathit{f}+{\mathit{f}}_{\mathit{y},\mathit{y}}{\mathit{f}}^2\right)(t_n,{\mathit{y}}_n)+\mathcal{O}\left(\delta t^4\right).\hfill \end{array}$$

The coefficients $a_1$ and $a_2$ satisfy $a_1+a_2=1$ and $2a_1{a}_2=1$. Thus, one can write ${\mathit{y}}_{n+1,2}=Re\left({\mathit{y}}_{n+1,2}\right)+iIm\left({\mathit{y}}_{n+1,2}\right)$. The real part of the approximated solution reads

$$\begin{array}{cc}\hfill Re\left({\mathit{y}}_{n+1,2}\right)& =Re\left({\mathit{y}}_n+\delta t{\mathit{y}}^{\prime }(t_n)+\frac{1}{2}\delta t^2{\mathit{y}}^{″}(t_n)+\frac{1}{4}\delta t^3{a}_1\left({\mathit{f}}_{t,t}+2{\mathit{f}}_{t,\mathit{y}}\mathit{f}+{\mathit{f}}_{\mathit{y},\mathit{y}}{\mathit{f}}^2\right)(t_n,{\mathit{y}}_n)\right.\hfill \\ \hfill & +\mathcal{O}\left(\delta t^4\right))\hfill \\ \hfill & =Re\left({\mathit{y}}_n+\delta t{\mathit{y}}^{\prime }(t_n)+\frac{1}{2}\delta t^2{\mathit{y}}^{″}(t_n)\right)\hfill \\ \hfill & +Re\left(a_1\right)\frac{\delta t^3}{4}\left({\mathit{f}}_{t,t}+2{\mathit{f}}_{t,\mathit{y}}\mathit{f}+{\mathit{f}}_{\mathit{y},\mathit{y}}{\mathit{f}}^2\right)(t_n,{\mathit{y}}_n)+Re\left(\mathcal{O}\left(\delta t^4\right)\right)\hfill \\ \hfill & ={\mathit{y}}_n+\delta t{\mathit{y}}^{\prime }(t_n)+\frac{\delta t^2}{2}{\mathit{y}}^{″}(t_n)+\frac{\delta t^3}{8}\left({\mathit{f}}_{t,t}+2{\mathit{f}}_{t,\mathit{y}}\mathit{f}+{\mathit{f}}_{\mathit{y},\mathit{y}}{\mathit{f}}^2\right)(t_n,{\mathit{y}}_n)\hfill \\ \hfill & +Re\left(\mathcal{O}\left(\delta t^4\right)\right).\hfill \end{array}$$

By using the Taylor expansion of the exact solution $\mathit{y}(t_{n+1})$ at the vicinity of $t_n$ and the expressions of the higher derivatives of the solution, we evaluate the error between the exact solution and $Re\left({\mathit{y}}_{n+1,2}\right)$ as follows

$$\begin{array}{cc}\hfill e(t_n+\delta t)& =\mathit{y}(t_{n+1})-Re\left({\mathit{y}}_{n+1,2}\right)\hfill \\ \hfill & =\frac{\delta t^3}{24}{\mathit{y}}^{(3)}(t_n)+\frac{\delta t^3}{8}\left({\mathit{f}}_t{\mathit{f}}_{\mathit{y}}+{\mathit{f}}_{\mathit{y}}^2\mathit{f}\right)(t_n,{\mathit{y}}_n)+Re\left(\mathcal{O}\left(\delta t^4\right)\right).\hfill \end{array}$$

(8)

As a consequence, the composed scheme is indeed convergent and of second-order, provided $\mathit{f}$ and its derivatives are bounded. Furthermore, the imaginary part of ${\mathit{y}}_{n+1,2}$ is expressed as follows:

$$\begin{array}{cc}\hfill Im\left({\mathit{y}}_{n+1,2}\right)& =Im\left(a_1\right)\frac{\delta t^3}{4}\left({\mathit{f}}_{t,t}+2{\mathit{f}}_{t,\mathit{y}}\mathit{f}+{\mathit{f}}_{\mathit{y},\mathit{y}}{\mathit{f}}^2\right)(t_n,{\mathit{y}}_n)+Im\left(\mathcal{O}\left(\delta t^4\right)\right)\hfill \\ \hfill & =-\frac{\delta t^3}{8}{\mathit{y}}^{(3)}(t_n)+\frac{\delta t^3}{8}\left({\mathit{f}}_t{\mathit{f}}_{\mathit{y}}+{\mathit{f}}_{\mathit{y}}^2\mathit{f}\right)(t_n,{\mathit{y}}_n)+Im\left(\mathcal{O}\left(\delta t^4\right)\right).\hfill \end{array}$$

By taking the difference between the error in (9) and the imaginary part of ${\mathit{y}}_{n+1,2}$, we end up with

$$e(t_n+\delta t)-Im\left({\mathit{y}}_{n+1,2}\right)=\frac{\delta t^3}{6}{\mathit{y}}^{(3)}(t_n)+\mathcal{O}\left(\delta t^4\right),$$

As a matter of fact, the imaginary part denoted in the sequel by ${\zeta }_{n+1}=Im\left({\mathit{y}}_{n+1,2}\right)$, is an approximation of the error up to third order and represents an error estimate of the approximated solution ${\mathit{y}}_{n+1}=Re\left({\mathit{y}}_{n+1,2}\right)$.    □

3.2. Time-Stepping Strategy

The vesicle can undergo more or less rapid dynamics or deformations over time. Although using a uniformly large time step can produce unphysical solutions and numerical instabilities, using larger time steps for some slow dynamics (e.g., around the horizontal position during tumbling) improves the computational efficiency. In [33], a manual adjustment is introduced to improve the quadratic convergence of the Newton’s algorithm. In this work, we adopt an adaptive time stepping scheme where the time step size is flexibly controlled such that the time integration error, estimated by ${\zeta }_n$ at $t_n$, does not exceed a prescribed threshold value. Let ${\delta t}_{max}$ and ${\delta t}_{min}$ be the upper and lower bounds of $\delta t$, generally chosen as ${10}^{-1}$ and 10 ${}^{-3}$ for the vesicle problem. Let C be a predefined parameter introduced to calibrate the level of adaptability, while ${10}^{-\alpha }$ with $\alpha >0$ is the desired precision and p represents the order of the method. The variable time step ${\delta t}_{n+1}$ reads:

$${\delta t}_{n+1}={\delta t}_n{\left(\frac{{10}^{-k}}{C\parallel {\zeta }_{n+1}{\parallel }_{L^2(\mathsf{\Lambda })}}\right)}^{1/p}.$$

Following [46], we also perform a comparison with a known criterion based on the the change rate of two consecutive numerical solutions ${\partial }_{\delta t}{\phi }_n$. For a predetermined constant $C={10}^3$, the adapted step size is given by:

$${\delta t}_{n+1}=max\left\{{\delta t}_{min},\frac{{\delta t}_{max}}{\sqrt{1+C\parallel {\partial }_{\delta t}{\phi }_n{\parallel }_{L^2(\mathsf{\Lambda })}^2}}\right\}.$$

3.3. Overall Algorithm

The membrane and fluid problems are first solved iteratively, while a fixed-point loop allows a strongly coupled approach for better numerical stability. Indeed, the level set problem is first solved using the velocity in the previous iteration; the fluid problem is subsequently solved using the updated level set function in the membrane force calculation. Fixed point convergence is reached when the absolute error between two subsequent iterates is less than a prescribed tolerance, set $\mathtt{tol}={10}^{-6}$ in the calculations. The overall coupling approach is detailed in Algorithm 1.

Algorithm 1 Fluid/vesicle coupling algorithm

1: Set fixed-point tolerance $\mathtt{tol}$ and time step $\delta t_0$ 2: Set initial conditions ($n=0$): ${\mathbf{u}}_0$ and ${\phi }_0$ are known values 3: for $n=0,\cdots ,n_{\max }$ do 4:    Update time $t_{n+1}=t_n+\delta t_n$ 5:    Set, from known values, ${\mathbf{u}}_{n+1}^0={\mathbf{u}}_n$ and ${\phi }_{n+1}^0={\phi }_n$ 6:    for $k=0,\cdots ,k_{\max }$ do 7:      Compute ${\phi }_{n+1}^{k+1}$ from (7d), advected with velocity ${\mathit{u}}_{n+1}^k$ 8:      Compute the mixed variable $\mathsf{\Psi }$ (7c) and membrane force in (7a) with ${\phi }_{n+1}^{k+1}$ 9:      Compute ${\mathbf{u}}_{n+1}^{k+1}$ and $p_{n+1}^{k+1}$ by solving the linear generalised Stokes subproblem (7a) and (7b) 10:      Compute absolute error ${\mathtt{err}}_{k+1}:={∥{\mathcal{H}}_{\epsilon }\left({\phi }_{n+1}^{k+1}\right)-{\mathcal{H}}_{\epsilon }\left({\phi }_{n+1}^k\right)∥}_{L_2(\mathsf{\Omega })}$ 11:      if ${\mathtt{err}}_{k+1}<\mathtt{tol}$ then 12:         Update solution ${\mathbf{u}}_{n+1}={\mathbf{u}}_{n+1}^{k+1}$, $p_{n+1}=p_{n+1}^{k+1}$ and ${\phi }_{n+1}={\phi }_{n+1}^{k+1}$ 13:         break 14:      end if 15:    end for 16:    if $t_{n+1}⩾T$ then 17:      break 18:    end if 19:    Compute the adapted time step $\delta t_{n+1}$ with ${\phi }_{n+1}$ 20: end for

4. Numerical Results and Discussion

We provide several validation tests to illustrate the accuracy in terms of spatial and temporal convergence and the ability of the model to predict experimentally observed outcomes, including vesicle regimes under simple shear flow in purely Newtonian flow.

The simulations have been run using the C++ library for scientific computing Rheolef [47], with particular emphasis on finite elements and parallel computing. The calculations were performed on non-regular meshes generated using Gmsh (GMSH– http://www.geuz.org/gmsh (accessed on 28 March 2023)). Figures are generated using Paraview (Paraview–http://www.paraview.org (accessed on 28 March 2023)) and Gnuplot (Gnuplot–http://www.gnuplot.info (accessed on 28 March 2023)).

4.1. Verification of the Method–Time Convergence

In this test case, we use numerical simulations to verify the above analytical results, mainly in terms of convergence behavior. We first test the accuracy of the time integration scheme and the composition method by solving the following initial value problem

$$y^{\prime }+tan(t-2)y=\frac{2a(t-2)cos(t-2)}{{(a(t-1)(t-3))}^2+1},y(0.5)=y_0.$$

The exact solution

$$y(t)=cos(t-2)arctan\left(a(t-1)(t-3)\right)$$

has stiff variations at $t=1$ and $t=3$ which are calibrated by the parameter a. The initial condition $y_0$ is given by the exact solution evaluated at $t=0.5$. We solve the differential equation numerically over the time interval $[0.5,6]$ using different numerical schemes for a relatively large stiffness coefficient $a=50$.

We consider constant time step sizes and perform a convergence study for the first-order basic BE, the double composition of BE, the second-order Cranck Nicolson (simply denoted by CN) versus successively refined step sizes $\delta t$. The log-log plot of errors with respect to the exact solution in the norm $l_{\infty }$ is reported in Figure 2 as a function of $\delta t$. Results corroborate the theoretical time rates of convergence, showing first-order and second-order convergences for BE and CN, respectively, while a higher quadratic convergence is achieved by dual composition of BE scheme. Moreover, a high order cubic convergence is obtained by double composition of the composed BE scheme. Consequently, the error bounds are of optimal order in time.

We now proceed to a quantitative comparison between the time-stepping strategies. The initial value problem is solved using adapted time steps generated using both the error estimate ${\zeta }_n$ and the rate of change criterion. The calculations are also performed with a constant time step $\delta t^{*}$ so as to perform the same number of iterations as in the case of the adaptation of $\delta t$ with the criterion ${\zeta }_n$. Results are reported in Figure 3a–c.

Figure 3c highlights the need to use an adapted time step according to the variations of the solution in order to obtain an accurate evaluation in areas of stiff variations. Visually, one perceives significant discrepencies compared to the exact solution, which is not the case for calculations with the adapted time steps.

We give on Figure 3a the evolution over time of the adapted time steps ${\delta t}_n$ using for the adaptation the two criteria above: the error estimate ${\zeta }_n$ and rate of change $|{\partial }_{\delta t}(y_n)|$. Although there is overall agreement in terms of variations of the adapted time steps, the ${\zeta }_n$ criterion makes it possible to capture the stiff variations with an appropriate level of temporal adaptability compared to the rate of change. Indeed, we plot on Figure 3c the error with respect to the exact solution as a function of time, showing clearly better accuracy when using ${\zeta }_n$ compared to the criterion $|{\partial }_{\delta t}(y_n)|$, which provides better results compared to the compuations with a constant time step $\delta t^{*}$.

4.2. Level Set Example–Spatial Convergence

To investigate spatial accuracy, we consider the Zalesak test of a solid slotted disk in a periodic rotation. This test is commonly used to assess the performance of interface tracking methods [48]. Consider a slotted circle of radius $0.15$ initially centered at $(0.5,0.75)$ and having a slot of depth $1/4$ and width $1/10$. The slotted circle is advected by a rotational field given by $\mathit{u}={\left({\displaystyle \frac{1}{2}(0.5-y),\frac{1}{2}(x-0.5)}\right)}^T$. The computational domain is the square ${[0,1]}^2$. Successive snapshots showing the rotating slotted disk are provided in Figure 4.

We perform a spatial convergence analysis by calculating the error in the L2 norm of the Heaviside of $\phi$ with respect to the exact solution after one period of rotation. Successively refined meshes with mesh size h are considered. Figure 5 plots the evolution of the logarithm of errors as a function of the logarithm of h, for continuous piecewise polynomial approximations of degrees $k=1$, $k=2$ and $k=3$. A nonlinear fitting is done. Convergence slopes show that the error behaves approximately as $\mathcal{O}\left(h^{k+1}\right)$, where k is the Lagrange polynomial’s degree.

4.3. Vesicle in Linear Shear Flow

Consider a vesicle with a reduced area $v=0.77$ in a simple linear shear flow and set the physical parameters to $Re={10}^{-3}$, $\mathrm{Bn}={10}^3$ and $\varsigma =0.479$. In light of the experimental observations and analysis [49], RBCs can undergo mainly two different motions, referred to as tank-treading (TT) and tumbling (TB) motions.

We first set a viscosity ratio $\overline{\mu }=1$ and report on Figure 6 successive snapshots of the membrane dynamics in flow. Observe that the membrane $\mathsf{\Gamma }$ follows a steady state TT regime by fixing its main axis at a given equilibrium inclination angle with respect to the flow. Isocontour lines are provided in Figure 6 showing the velocity field becomes tangential to the membrane in the steady state. Indeed, each material point on the membrane continues to rotate tangentially to the flow without the vesicle changing shape.

We then increase the ratio to $\overline{\mu }=10$, which leads to a change of regime. The membrane loses its fluidity and undergoes a solid-like cyclic rotational motion around the interior, as shown in the snapshots in Figure 7. Indeed, the TT becomes unfavorable to the benefit of the TB when the viscosity ratio is beyond a threshold value for which the twisting force due to shear can no longer be converted into TT torque of the membrane.

Figure 8 reports the TB inclination angle ${\theta }_n$ and adapted time steps versus the time for few rotation periods. Higher values of ${\delta t}_n$ are successfully generated when the vesicle is close to a horizontal tilt characterized by slow rotational dynamics. However, ${\delta t}_n$ are almost ten times smaller when the membrane is in vertical position undergoing a strong action of the orthogonal shear flow which results in a faster rotation.

Thereafter, we proceed to a qualitative validation of our model compared to some experimental and numerical results available in the published literature. We first focus on the TT-TB transition in the purely Newtonian case, where some simulation results are available. We report the critical viscosity contrast ${\overline{\mu }}^{\star }$ necessary to observe a change of regime for several values of the reduced area v in the interval $[0.7,1]$. The phase diagram is reported in Figure 9, showing good quantitative agreement with the Keller and Skalak’s 2D predictive theory [50], referred to as KS, the experimental results of Kantsler & Steinberg in [51] and the numerical calculations with the level set method in a finite differences framework in [26]. Note that the KS predictive theory is obtained by assuming a fixed ellipse shape for the vesicle, while only flow orientation is permitted. This could explain the small difference between the different results and those predicted by the KS work.

Finally, we consider a vesicle following the tank-treading motion and focus on some effects of the non-Newtonian fluid model on the TT regime. For TT vesicles at equilibrium, we report in Figure 10 the change in the equilibrium inclination angle ${\theta }^{\star }$ with respect to the reduced area parameter v. For $\varsigma =1$ corresponding to a Newtonian fluid, only a qualitative agreement is observed with the experimental results in [51] and the numerical calculations in [12], with angles slightly shifted upwards. While the proposed penalty method features a smaller system of equations to solve, the small difference observed with respect to the numerical results obtained in [12,33] may be due to relaxation of the inextensibility instead of using an exact Lagrange multiplier [18,33] imposing this constraint.

We now study the equilibrium inclination angle of tank-treading vesicles immersed in a non-Newtonian Carreau fluid. Simulation results are also reported in Figure 10, illustrating lower values of the equilibrium angle ${\theta }^{\star }$ are obtained, for various reduced area parameters. Our interpretation is as follows. At the vesicle scale in small capillaries, the length scale and velocity are typically around $10\mathsf{\mu }$m and 1–100 $\mathsf{\mu }\mathrm{m}/\mathrm{s}$, respectively. The local shear rate is of the order of 0.1–10 $<100{\mathrm{s}}^{-1}$ (typical value characterizing the transition to a Newtonian behavior), resulting in an apparent viscosity closer to the infinite-shear viscosity. As a consequence, the fluid is much denser and tends to decrease the inclination angle at equilibrium.

5. Conclusions

We have presented a numerical method for the simulation of the dynamics of a single vesicle in non-Newtonian flow. The membrane is described by a level set function, while the inextensibility is relaxed by a penalty method. The penalty method is straightforward to implement using any Stokes solver generally available in finite element libraries and features a smaller system of equations to solve. The order of the temporal scheme is raised according to the double composition of the numerical flow associated with the first-order backward Euler integrator. Moreover, an adaptive time-stepping strategy is based on an error estimation and allows to improve the prediction of the membrane dynamics. The numerical results show that we are able to accurately reproduce the cellular regimes under linear shear flow as well as the regime change from tank-treading to tumbling in the Newtonian case. Last but not least, a change of regime due to a non-Newtonian rheology is highlighted.

The non-Newtonian effect will be further explored in a separate work, in the hope of prompting experiments to explore it in depth. This is part of a larger study of the dynamics of 3D biological membrane and RBCs with hyperelastic cytoskeleton in non-Newtonian flows [25,52]. This requires the development of a robust solver for fluid–structure interaction problems with a thin elastic structure [53,54]. Other extensions of the presented method are being currently explored. We are particularly focusing on the development of higher-order schemes by multi-step composition of various basic integrators. The construction of robust preconditioners is also being explored.