# A Finite Element Approach For Modeling Biomembranes In Incompressible Power-Law Flow

Aymen Laadhari, Ahmad Deeb

<figure><img src="image_1.png" /></figure>

## To cite this version:

Aymen Laadhari, Ahmad Deeb. A Finite Element Approach For Modeling Biomembranes In Incompressible Power-Law Flow. 2023. hal-04102960

HAL Id: hal-04102960

https://hal.science/hal-04102960

Preprint submitted on 22 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Finite Element Approach For Modeling Biomembranes In Incompressible Power-Law Flow

Aymen Laadhari<sup>a)</sup> and Ahmad Deeb

Department of Mathematics, College of Arts and Sciences, Khalifa University of Science and Technology, Abu Dhabi, United Arab Emirates.

<sup>a)</sup>Corresponding author: Aymen.Laadhari@ku.ac.ae

**Abstract.** We present a numerical method to model the dynamics of inextensible biomembranes in a quasi-Newtonian incompressible flow, which better describes hemorheology in the small vasculature. We consider a level set model for the fluid-membrane coupling, while the local inextensibility condition is relaxed by introducing a penalty term. The penalty method is straightforward to implement from any Navier-Stokes/level set solver and allows substantial computational savings over a mixed formulation. A standard Galerkin finite element framework is used with an arbitrarily high order polynomial approximation for better accuracy in computing the bending force. The PDE system is solved using a partitioned strongly coupled scheme based on Crank-Nicolson time integration. Numerical experiments are provided to validate and assess the main features of the method.

## INTRODUCTION

This paper is concerned with the numerical study of the time-dependent dynamics of biomembranes in a surrounding Newtonian and non-Newtonian flow. The coupled fluid-membrane problem is highly nonlinear and time consuming.

Blood is a very complex fluid. Its rheology at the macroscopic scale depends both on the individual dynamics of its embedded entities and their fluid-structure interactions at the microscopic level. Red blood cells, referred to as RBCs, represent its main cellular component; They are responsible for the supply of oxygen and the capture of carbon dioxide. In the laboratory, giant unilamellar vesicles (diameter ≈ 10µm) are biomimetic artificial liquid drops, used in vitro and in silico to study the RBCs. Understanding the dynamics of RBCs in flow remain a difficult problem in the field of computational physics and at the theoretical level as well, consequently leading to a growing interest in the past two decades. In the published literature, several works have covered the areas of experimental biology [1], theoretical biology [2], physics [3, 4, 5] and applied mathematics [6, 7].

From a mechanical continuum perspective, Canham [8], Helfrich [9] and Evans [10] independently introduced in the early 1970s a model to describe the mechanics of lipid bilayer membranes, where cellular deformations are driven by the principal curvatures. This results in a highly nonlinear membrane force with respect to shape, see a mathematical derivation for a generalized energy functional based on shape optimization in [11].

Different methods have been developed to study the dynamics of biomembranes in a Newtonian flow. We can distinguish the level set method [12, 13, 14], the phase field method [15], the immersed boundary method [16], the boundary integral method [17], the parametric finite elements [7], and the lattice Boltzmann method [18]. From a numerical point of view, iterative and fully explicit decoupling strategies for the membrane-fluid problem are the most used techniques [14, 19]. An explicit treatment of the bending force usually leads to numerical instability problems and severe time step limitations, depending on the local mesh size and bending stiffness. However, only few works devised semi-implicit [7] or fully implicit time integration schemes [13, 20]. Although stability is improved, a high computational burden is generally obtained with implicit strategies. Other interesting decoupling strategies can be found in [21, 22, 23, 24].

While blood flow behaves like Newtonian fluid in larger diameter arteries at high shear rates, it exhibits non-Newtonian behavior in small diameter arteries with low shear rates at the microscopic scale [25]. Non-Newtonian rheology is mainly due to polymerization and the underlying mechanisms leading to the activation and deactivation of platelets and the interactions between different microscopic entities. Blood viscosity tends to increase at low shear rates as RBCs aggregate into a roller shape. The Casson, Power-Law, and Quemada models are the most widely used generalised Newtonian rheologies for blood [26, 27]. To our knowledge, such models have not yet been studied for the current problem. In this work, we consider a quasi-Newtonian power law model to describe the hemorheology.

The aim of this paper is to study the dynamics of biomembranes in a complex non-Newtonian incompressible viscous flow. In order to keep a reasonable computational cost compared to a fully mixed formulation, we design a penalty method to account for the local inextensibility of the membrane. Various higher-order finite element approximations are used to better approximate the bending force. We present a set of numerical examples to validate and show the main features of the method.
<figure><img src="image_3.png" /><figcaption>FIGURE 1: Sketch of the membrane $\Gamma$ embedded into a computational domain $\Lambda$, while $\Omega$ is the inner region.</figcaption></figure>

# MATHEMATICAL SETTING

## Membrane model

The deformations of the membrane allow minimizing the Canham-Helfrich-Evans [8, 9] bending energy while preserving the local inextensibility of the membrane. Let $H$ be the mean curvature, corresponding to the sum of the principal curvatures on the membrane. In the two-dimensional case, the membrane minimizes the bending energy given by:

$$J(\Omega) = \frac{k_b}{2} \int_{\partial\Omega} (H(\Omega))^2 ds, \qquad (1)$$

where $k_b \approx 10^{-20}/10^{-19}kgm^2s^{-2}$ is the bending rigidity modulus. The energy is a variant of the Willmore energy [28]. Let $T$ be the final time of the experiment. For any time $t \in [0, T]$, $\Omega(t) \subset \mathbb{R}^d$, $d = 2, 3$, is the interior domain of the membrane $\Gamma(t) = \partial\Omega(t)$, assumed Lipschitz continuous. The membrane is embedded in the domain $\Lambda$ which is large enough so that $\Gamma(t) \cap \partial\Lambda = \emptyset$, see Fig. 1. Hereafter, the dependence of $\Omega$ and $\Gamma$ upon $t$ is dropped to alleviate notations.

For a membrane with fixed topology, the Gauss-Bonnet theorem [29] states that the energy term weighted by $k_g$ is constant and can be ignored. The spontaneous curvature helps describe the asymmetry of phospholipid bilayers at rest, e.g. when different chemical environments exist on either side of the membrane. We assume $H_0 = 0$. Let $n$ and $\nu$ be the outward unit normal vector $\Gamma(t)$ and on $\partial\Lambda$, respectively. We introduce the surface gradient $\nabla_s = (\mathbf{Id} - n \otimes n) \nabla$, surface divergence $\text{div}_s = \text{tr}(\nabla_s)$ and surface Laplacian $\Delta_s = \text{div}_s(\nabla_s)$, where $\mathbf{Id}$ is the identity tensor. The expression and derivation of the bending force using shape optimization tools can be found in [11].

Membrane deformations are subject to specific constraints. Fluid incompressibility is assumed, this is $\text{div} u = 0$ in $\Lambda$. In addition, RBCs are phospholipid bilayers with local membrane inextensibility. This corresponds to a zero surface divergence, i.e. $\text{div}_s u = 0$ over $\Gamma$, that helps preserve the local perimeter. Global perimeter conservation follows from Reynolds' lemma [13]. As a consequence, a saddle point formulation results in a membrane surface force that balances the jump in hydrodynamic stress tensor and appears in the right side of (3g).

## Level set description

The motion of the membrane is followed implicitly in a level set framework as the zero level set of a function $\varphi$. For $t \in ]0, T[$, $\varphi$ is initialized by a signed distance $\varphi_0$ to $\Gamma(0)$ and satisfies the transport equation (3a), with $u$ the advection vector and $\varphi = \varphi_b$ on the upstream boundary $\Sigma_- = \{x \in \partial\Lambda : u \cdot \nu(x) < 0\}$. Geometric quantities such as $n = \nabla\varphi/|\nabla\varphi|$, $H = \text{div}_s n$ and bending force are coded in terms of $\varphi$ and are then extended to the entire computational domain $\Lambda$. Over time, a redistancing problem is resolved to maintain the signed distance property lost
by advection [30]. Indeed, a too large or too small gradient of $\varphi$ close to $\Gamma$ deteriorates the precise computation of the surface terms. Let $\varepsilon$ be a regularization parameter. We introduce the regularized Heaviside $\mathcal{H}_\varepsilon$ and Dirac $\delta_\varepsilon$ functions:

$$ \mathcal{H}_\varepsilon(\varphi) = \begin{cases} 0, & \text{when } \varphi < -\varepsilon \\ \frac{1}{2}\left(1 + \frac{\varphi}{\varepsilon} + \frac{1}{\pi}\sin\left(\frac{\pi\varphi}{\varepsilon}\right)\right), & \text{when } |\varphi| \le \varepsilon, \\ 1, & \text{otherwise} \end{cases} \quad \text{and} \quad \delta_\varepsilon(\varphi) = \frac{\mathrm{d}\mathcal{H}_\varepsilon}{\mathrm{d}\varphi}(\varphi). $$

Given a function $\zeta$ defined on $\Gamma$ and its extension $\tilde{\zeta}$ to $\Lambda$, surface integrals are approximated as follows:

$$ \int_{\Gamma} \zeta(x) \, ds \approx \int_{\Lambda} |\nabla \varphi| \delta_{\varepsilon}(\varphi) \tilde{\zeta}(x) \, dx. $$

## Governing equations

We assume constant densities $\rho_i$ and $\rho_o$ inside and outside of the membrane, respectively. Let us introduce the fluid velocity $u$ and the pressure $p$ which represent a Lagrange multiplier corresponding to the incompressibility constraint on $\Lambda$. Analogously, a position-dependent surface tension $\lambda$ helps imposing the local inextensibility constraint on $\Gamma$. Let $\mathbf{D}(u) = (\nabla u + \nabla u^T)/2$ be the shear strain rate tensor, so the fluid Cauchy stress tensor is $\sigma = \mathbf{T} - p\mathbf{I}$ where $\mathbf{T}$ is the stress deviator. The normal stress jump $[\sigma n]_\perp = \sigma_+ n - \sigma_- n$ on $\Gamma$ describes the interactions of the membrane with the surrounding fluid [20], while the stress discontinuity is calibrated by (3f). For a simple shear flow, $u_b$ is the shear rate on $\Sigma_D \subset \partial\Lambda$, while natural boundary conditions are prescribed on $\Sigma_N \subset \partial\Lambda$.

We assume a quasi-Newtonian power-law model [26] where the nonlinear constitutive equation expresses the stress deviator with a power-law viscosity function as

$$ \mathbf{T} = 2\eta (|\mathbf{D}(u)|^2) \mathbf{D}(u), \quad \text{with } \eta(\gamma) = K\gamma^{(\nu-1)/2}, \quad \text{for all } \gamma \in \mathbb{R}, \qquad (2) $$

where $\nu > 0$ and $K$ are the power index and consistency index, respectively. According to [31], $\nu = 0.7755 < 1$ (i.e. a shear thinning fluid) and $K = 14.67 \times 10^{-3}$ Pa s for normal blood samples obtained using a multiple regression technique. The Newtonian case $\nu = 1$ corresponds to a linear stress-strain relationship that reduces the viscosity function $\eta(\gamma) = K$ to a constant. By analogy with the Newtonian case, $K = \mu_i$ and $K = \mu_o$ stand for the values of the consistency index in the intra- and extra-membrane domains, respectively.

We perform a dimensionless analysis. Let $U$ be the maximum velocity on $\Sigma_D$ and $D$ the diameter of a circle having the same membrane perimeter. We consider the dimensionless Reynolds number $Re = \rho_o U D \mu_o^{-1}$ which expresses the ratio between the inertial and viscous forces, and the capillary number $Ca = \mu_o D^2 U k_b^{-1}$ which compares the flow force to the bending resistance of the membrane. Furthermore, the parameter $\beta = \mu_i/\mu_o$ represents the ratio of consistency indices and corresponds to the viscosity ratio with respect to extracellular viscosity in the Newtonian case. The regularized dimensionless viscosity function is:

$$ \mu_{\varepsilon}(\varphi)|\mathbf{D}(u)|^{\nu-1} = (\mathcal{H}_{\varepsilon}(\varphi) + \beta(1 - \mathcal{H}_{\varepsilon}(\varphi)))|\mathbf{D}(u)|^{\nu-1}. $$

Following [20], we choose $\rho_i = \rho_o$. Let $\sigma_\varepsilon$ stand for the regularized Cauchy stress tensor. The dimensionless reduced area $\Xi_{2d} = 4\pi|\Omega|/|\Gamma|^2 \in ]0, 1]$ compares the area of the interior domain to that of a circle with the same perimeter. The dimensionless coupled problem writes: find $\varphi$, $u$, $p$ and $\lambda$ such that

$$ \partial_t \varphi + \mathbf{u} \cdot \nabla \varphi = 0 \quad \text{in } ]0, T[ \times \Lambda \qquad (3a) $$

$$ \mathrm{Re} (\partial_t u + \mathbf{u} \cdot \nabla u) - \mathbf{div}(\sigma_{\varepsilon}(D(u), p, \varphi)) = 0 \quad \text{in } ]0, T[ \times (\Lambda \setminus \partial \Omega) \qquad (3b) $$

$$ \mathrm{div} \, \mathbf{u} = 0 \quad \text{in } ]0, T[ \times \Lambda \qquad (3c) $$

$$ \mathrm{div}_s \mathbf{u} = 0 \quad \text{on } ]0, T[ \times \partial \Omega \qquad (3d) $$

$$ [\mathbf{u}]_{-}^{+} = 0 \quad \text{on } ]0, T[ \times \partial \Omega \qquad (3e) $$

$$ [\sigma_{\varepsilon} \mathbf{n}]_{-}^{+} = \nabla_s \lambda - \lambda H \mathbf{n} + (2Ca)^{-1} (2\Delta_s H + H^3) \mathbf{n} \quad \text{on } ]0, T[ \times \partial \Omega \qquad (3f) $$

$$ \varphi = \varphi_b \quad \text{on } ]0, T[ \times \Sigma_{-} \qquad (3g) $$

$$ \mathbf{u} = \mathbf{u}_b \quad \text{on } ]0, T[ \times \Sigma_D \qquad (3h) $$
$$ \sigma \cdot \nu = 0 \quad \text{on } ]0, T[ \times \Sigma_N \qquad (3i) $$

$$ \varphi(0) = \varphi_0 \quad \text{in } \Lambda \qquad (3j) $$

$$ \mathbf{u}(0) = \mathbf{u}_0 \quad \text{in } \Lambda. \qquad (3k) $$

Let $\varepsilon_\lambda = 10^{-8}$ be the penalty parameter. To make the method straightforward to implement from any Level Set / Navier-Stokes solver and considerably reduce the size of the linear system to be solved, the inextensibility constraint is relaxed by introducing a penalty term. Indeed, the corresponding minimization problem should be approximated by another minimization problem by penalizing the local inextensibility constraint for the velocity (3d). See analogous penalty method for other applications in [32].

To overcome instability problems when solving the level set equation using the standard Galerkin method, there are a variety of stabilization methods such as the streamline diffusion method, the subgrid viscosity method and the Streamline Upwind Petrov-Galerkin (SUPG) method used in this work. The latter introduces a stabilization term by adding a diffusion in the streamline direction.

We introduce the functional spaces of admissible velocity $\mathbf{u}$, pressure $p$ and level set $\varphi$:

$$ \mathcal{V}(\mathbf{u}_b) = \{v \in (H^1(\Lambda))^d : v = \mathbf{u}_b, \text{ on } \Sigma_D\}, \quad \mathcal{Q} = \{q \in L^2(\Lambda) : \int_{\Omega} q = 0\}, $$

$$ \mathcal{X}(\varphi_b) = \{\psi \in W^{1,\infty}(\Lambda) \cap H^1(\Lambda) : \psi = \varphi_b, \text{ on } \Sigma_-\}. $$

To reduce a derivation order of $\varphi$ when evaluating the bending strength, we use the Green formula on a closed surface. See e.g. [13]. Testing with appropriate test functions and integrating (3b) over $\Omega$ and $\Lambda \setminus \bar{\Omega}$ separately, the variational problem writes:

Find $\mathbf{u} \in \mathcal{C}^0([0, T[, L^2(\Lambda)^d) \cap L^2([0, T[, \mathcal{V}(\mathbf{u}_b))$, $p \in L^2([0, T[, \mathcal{Q})$, and $\varphi \in \mathcal{C}^0([0, T[, L^2(\Lambda)^d) \cap L^2([0, T[, \mathcal{X}(\varphi_b))$ such that

$$ \mathrm{Re} \int_{\Lambda} \left( \frac{\partial \mathbf{u}}{dt} + \mathbf{u} \cdot \nabla \mathbf{u} \right) \cdot \mathbf{v} + \int_{\Lambda} 2\mu_{\varepsilon}(\varphi) |\mathbf{D}(\mathbf{u})|^{v-1} \mathbf{D}(\mathbf{u}) : \mathbf{D}(\mathbf{v}) + \frac{1}{\varepsilon_{\lambda}} \int_{\Lambda} \mathrm{div}_s(\mathbf{u}) \mathrm{div}_s(\mathbf{v}) |\nabla \varphi| \delta_{\varepsilon}(\varphi) \\ - \int_{\Lambda} p \, \mathrm{div} \, \mathbf{v} + \frac{1}{2\mathrm{Ca}} \int_{\Lambda} \delta_{\varepsilon}(\varphi) |\nabla \varphi| (2\nabla_s H \cdot \nabla_s (n \cdot \mathbf{v}) - H^3 n \cdot \mathbf{v}) = \int_{\Sigma_N} \sigma \mathbf{v} \cdot \mathbf{v}, \quad \forall \mathbf{v} \in \mathcal{V}(0), \quad (4a) $$

$$ \int_{\Lambda} q \, \mathrm{div} \, \mathbf{u} = 0, \qquad \forall q \in \mathcal{Q}, \qquad (4b) $$

$$ \int_{\Lambda} \frac{\partial \varphi}{\partial t} \psi + \int_{\Lambda} (\mathbf{u} \cdot \nabla \varphi) \psi + \int_{\Lambda} \xi(\tau; \varphi, \psi) = 0, \qquad \forall \psi \in \mathcal{X}(0). \quad (4c) $$

Here, $\xi(\tau; \varphi, \psi)$ stands for the SUPG stabilisation term and $\tau$ is a stabilization parameter defined element wise to control the amount of diffusion.

## NUMERICAL APPROACH

The interval $[0, T]$ is divided into $N$ sub-intervals $[t^n, t^{n+1})$ with $0 \le n \le N-1$ of constant step $\Delta t$. For $n > 0$, $\mathbf{u}^n$, $p^n$ and $\varphi^n$ are computed by induction to approximate $\mathbf{u}$, $\varphi$ and $p$ at $t^n$. We use the Crank-Nicolson scheme for the time discretization of (3a) and (3b) without the need to bootstrap the initial conditions. The choice of this scheme was for its simplicity to implement and being a second order one-step integrator. The discretized (3a) writes

$$ \varphi^{n+1} = \varphi^n + \frac{\Delta t}{2} (\mathbf{u} \cdot \nabla \varphi^{n+1} + \mathbf{u} \cdot \nabla \varphi^n) \quad \text{in } \Lambda. $$

For the spatial discretization, we consider a partition $\mathcal{T}_h$ of $\Lambda$ consisting of geometrically conformal open simplicial elements $K$. We define the mesh size as the diameter of the largest mesh element $h = \max h_K$ with $K \in \mathcal{T}_h$.

We consider a Taylor-Hood finite element approximation for $\mathbf{u}$ and $p$. After using a surface Green's transformation, the evaluation of the Canham-Helfrich-Evans force requires a third-order derivative in $\varphi$ which induces numerical oscillations when using lower-order polynomial approximations. To avoid introducing additional mixed variables and additional equations as in [13], higher degree polynomials are considered for the discretization of $\varphi$ because the bending force requires its fourth order derivatives. For the SUPG method, the streamline diffusion parameter is
chosen numerically proportional to the local mesh size, this is $\tau_K = Ch_K/\max\{|u|_{0,\infty,K}, \text{tol}/h_K\}$, where C is a scaling constant and tol/h$_K$ helps to avoid division by zero. To overcome the instability problems induced by an explicit decoupling, we consider a partitioned implicit strategy based on a fixed point algorithm, as detailed in Alg. 1.

<table><caption>**Algorithm 1** Fluid-membrane coupling</caption><tr><td>1:</td><td>n = 0: let φ<sup>0</sup> and u<sup>0</sup> being given</td></tr><tr><td>2:</td><td>for n = 0, ..., N − 1 do</td></tr><tr><td>3:</td><td>Initialize u<sup>n+1,0</sup> = u<sup>n</sup>, φ<sup>n+1,0</sup> = φ<sup>n</sup></td></tr><tr><td>4:</td><td>while e<sup>k</sup> &lt; 10<sup>−6</sup> do</td></tr><tr><td>5:</td><td>Compute φ<sup>n+1,k+1</sup> using u<sup>n+1,k</sup></td></tr><tr><td>6:</td><td>Compute u<sup>n+1,k+1</sup>, p<sup>n+1,k+1</sup> using φ<sup>n+1,k+1</sup></td></tr><tr><td>7:</td><td>Compute the error e<sup>k</sup> = |u<sup>n+1,k+1</sup> − u<sup>n+1,k</sup>|<sub>1,2,Λ</sub>/|u<sup>n+1,k</sup>|<sub>0,2,Λ</sub> + |φ<sup>n+1,k+1</sup> − φ<sup>n+1,k</sup>|<sub>0,2,Λ</sub>/|φ<sup>n+1,k</sup>|<sub>0,2,Λ</sub></td></tr><tr><td>8:</td><td>end while</td></tr><tr><td>9:</td><td>Update u<sup>n+1</sup> = u<sup>n+1,k+1</sup>, φ<sup>n+1</sup> = φ<sup>n+1,k+1</sup></td></tr><tr><td>10:</td><td>end for</td></tr></table>

## NUMERICAL EXAMPLES

<figure><img src="image_4.png" /><figcaption>FIGURE 2: Reversible vortex. (Left) Snapshots showing the interface deformations at $t \in \{0, 0.25, 0.57, 0.75, 0.875, 1\}$ with $h = 0.01$. (Right) Spatial convergence in $L^2$ norm for high-order finite element approximations.</figcaption></figure>

### Example 1: Reversible Vortex - Grid convergence.

Simulations were performed using FEniCSx [33]. To evaluate the capability of the level set solver for high-order finite elements, necessary afterwards for an accurate assessment of highly nonlinear bending force, we consider a reversible vortex test case featuring large deformations of the interface. The computational domain is $\Lambda = [0, 1]^2$. A circular interface of radius $R = 0.15$ initially centered at $(0.7, 0.7)$ is stretched into thin filaments which are coiled like a starfish by a vortex flow field. The deformations are periodic and the stretching of the membrane unravels before the interface regains its circular shape after a period at $t = T$. The maximal deformation $\psi$ occur at $t = T/2$, with $\psi = 3$ and $T = 1$ in numerical computations. Similar 2D and 3D test cases are widely used to test interface tracking methods. We follow LeVeque's test [34] (Example 9.5) and consider a velocity field at $x = (x, y)^T \in \Lambda$ given by

$$u(t, x) = (-2 \sin(\psi \pi x)^2 \sin(\psi \pi y) \cos(\psi \pi y) \cos(\pi t/T), 2 \sin(\psi \pi y)^2 \sin(\psi \pi x) \cos(\psi \pi x) \cos(\pi t/T))^T.$$
The spatial accuracy of the finite element numerical approximations is studied by computing the errors in $L^2(\Lambda)$ norm on successively refined meshes with respect to an exact reference solution $\pi_h \varphi$ at $t = T$, where $\pi_h$ represents the Lagrange interpolation operator. Errors are calculated after one stretching period. For $k$ the degree of the polynomial approximation, the time step $\Delta t = h^k$ is chosen small enough not to significantly influence the overall accuracy. Fig. 2 reports the convergence of calculated errors with respect to the mesh size for several polynomial finite element approximations. Convergence rates are also displayed, showing for instance an almost second-order accuracy for $k = 1$ and fifth-order accuracy for $k = 4$.

<figure><img src="image_5.png" /><figcaption>FIGURE 3: TT regime: Change in $\theta^*/\pi$ with respect to $\Xi_{2d}$ for a viscosity ratio $\beta = 1$. Comparisons with results from [19] and [20] ($\text{Re} = 10^{-3}$, $\text{Ca} = 100$), [35] ($\text{Ca} = 9$) and [36] ($\text{Ca} = 10$).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4: TT regime: Change in $\theta^*/\pi$ with respect to $\Xi_{2d}$ for $\beta = 2.7$. Comparison of non-Newtonian model with results in [20, 35] and measurements in [37].</figcaption></figure>

## Example 2: Dynamics of the biomembrane in Newtonian and quasi-Newtonian flows.

We first proceed to a quantitative validation with some experimental and numerical results available in the literature in the case of a purely Newtonian flow. We set $v = 1$, a viscosity contrast $\beta = 1$, $Ca = 10^2$ and $\text{Re} = 9 \times 10^{-3}$. More details on the physiological values of the Reynolds number at the level of RBCs are available in [19]. The membrane follows a tank-treading type movement, called TT, where it reaches a steady state characterized by a fixed angle of inclination; The surrounding fluid continues its rotation tangentially to the membrane. We consider different values of the reduced areas $\Xi_{2d} \in [0.6, 1]$, and calculate the angle of inclination at equilibrium $\theta^*$. Fig. 3 and Fig. 4 plot the change in $\theta^*/\pi$ against $\Xi_{2d}$ in both Newtonian and quasi-Newtonian cases for different values of the viscosity ratio $\beta$. The results are compared with those of Kraus et al. [36], Zhao et al. [35], Salac et al. [19] and Laadhari et al. [20], showing good overall consistency. However, note that the values obtained with $v = 0.7755$ fit slightly better compared to those of the Newtonian model which are a little higher than the other curves. This cannot be confirmed at all, given the setting of different dimensionless values such as $\text{Re}$ and $\text{Ca}$ in the different experiments. An in-depth study is in progress and will be the subject of a forthcoming work.
<figure><img src="image_7.png" /><figcaption>FIGURE 5: TT regime for $\Xi_{2d} = 0.68$, $\beta = 1$, $Ca = 4 \times 10^4$ and $Re = 9 \times 10^{-3}$ at $t \in \{0, 0.125, 0.25, 0.5, 1, 2\}$.</figcaption></figure>

Simulations are now performed using a different ratio $\beta^* = 2.7$ in the non-Newtonian case. We calculate the angle of inclination $\theta^*/\pi$ and compare with some numerical [35] and experimental [37] results available only for larger reduced areas. Fig. 3(right) shows close but slightly higher equilibrium angles when the shape of the membrane becomes close to a circle. The deviations can be mainly due to the non-Newtonian model, but also to the different values of the confinement levels and the boundary conditions used in the different works.

According to an experimental systematic study on individual individual red cells in a simple shear flow, a change in dynamics occurs when the viscosity ratio exceeds a critical value depending on the reduced area [38]. This is the tumbling regime, noted by TB, which is characterized by the periodic rotation of the membrane around its axis. A well-known empirical model was developed by Keller and Skalak [4]. This dynamics was obtained in the simulations with the non-Newtonian model, see Fig. 5 and Fig. 6 for the snapshots of the TT and TB dynamics obtained with the same set of parameters but with $\beta = 1$ and $\beta = 10$, respectively.

<figure><img src="image_8.png" /><figcaption>FIGURE 6: Snapshots showing a membrane in TB regime for $\Xi_{2d} = 0.68$, $\beta = 10$, $Ca = 4 \times 10^4$ and $Re = 9 \times 10^{-3}$, at times $t \in \{0, 0.13, 0.25, 0.5, 1.25, 2, 2.25, 2.38, 3, 4, 4.5, 4.7\}$, respectively.</figcaption></figure>

CONCLUSION

We have presented in this paper a relatively simple method for simulating the dynamics of an individual red blood cell, or inextensible biological membrane in general, in a surrounding incompressible non-Newtonian flow that better describes the hemorheology in small capillaries.

We validated our framework using high-order finite element approximations in the case of a membrane in a simple shear flow. Simulations have shown that the method is capable of capturing the basic cellular dynamics, namely the well-known tank treading and tumbling motions. This is part of a larger ongoing work to explore the dynamics of red blood cells in small capillaries, while accounting for cell elasticity [39] in non-Newtonian surrounding flow.

ACKNOWLEDGMENTS

The authors acknowledge financial support from KUST through the grant FSU-2021-027.
REFERENCES

1. H. Song and Y. Hashash, "Characterization of stress-strain behaviour of red blood cells (RBCs), part II: response of malaria-infected RBCs," *Inverse Probl. Sci. Eng.* **23**, 413–424 (2015).

2. S. A. Safran, *Statistical Thermodynamics of Surfaces, Interfaces and Membranes*, Frontier in Physics Vol 90 (Addison-Wesley Publishing Company, Reading, Massachusetts, 1994).

3. B. Kaoui, G. H. Ristow, I. Cantat, C. Misbah, and W. Zimmermann, "Lateral migration of a two-dimensional vesicle in unbounded poiseuille flow," *Phys. Rev. E* **77**, 021903 (2008).

4. S.R.Keller and R. Skalak, "Motion of a tank-treading ellipsoidal particle in a shear flow," *J. Fluid Mech.* **120**, 27–47 (1982).

5. W. Choi, J. Yi, and Y. Kim, "Fluctuations of red blood cell membranes: The role of the cytoskeleton," *Phys. Rev. E* **92**, 012717 (2015).

6. G. Dziuk, "Computational parametric Willmore flow," *Numer. Math.* **111**, 55–80 (2008).

7. J. W. Barrett, H. Garcke, and R. Nurnberg, "Numerical computations of the dynamics of fluidic membranes and vesicles," *Phys. Rev. E* **92**, 052704 (2015).

8. P. Canham, "The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell," *J. Theor. Biol.* **26**, 61–81 (1970).

9. W.Helfrich, "Elastic properties of lipid bilayers: theory and possible experiments," *Z. Naturforsch. C* **28**, 693–703 (1973).

10. E. A. Evans, "Bending resistance and chemically induced moments in membrane bilayers," *Biophys. J.* **14**, 923–931 (1974).

11. A. Laadhari, C. Misbah, and P. Saramito, "On the equilibrium equation for a generalized biological membrane energy by using a shape optimization approach," *Physica D* **239**, 1567–1572 (2010).

12. G.-H. Cottet, E. Maitre, and T. Milcent, "Eulerian formulation and Level-Set models for incompressible fluid-structure interaction," *Math. Model. Numer. Anal.* **42**, 471–492 (2008).

13. A. Laadhari, P. Saramito, and C. Misbah, "Computing the dynamics of biomembranes by combining conservative level set and adaptive finite element methods," *J. Comput. Phys.* **263**, 328–352 (2014).

14. V. Doyeux, Y. Guyot, V. Chabannes, C. Prud'homme, and M. Ismail, "Simulation of two-fluid flows using a finite element/level set method. Application to bubbles and vesicle dynamics," *J. Comput. Appl. Math.* **246**, 251–59 (2013).

15. Q. Du, C. Liu, and X. Wang, "A phase field approach in the numerical study of the elastic bending energy for vesicle membranes," *J. Comput. Phys.* **198**, 450–468 (2004).

16. W.-F. Hu, K. Y. Kim, and M.-C. Lai, "An immersed boundary method for simulating the dynamics of three-dimensional axisymmetric vesicles in Navier-Stokes flows," *J. Comput. Phys.* **257**, Part A, 670–686 (2014).

17. A. Rahimian, S. K. Veerapaneni, and G. Biros, "Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method," *J. Comput. Phys.* **229**, 6466–6484 (2010).

18. B. Kaoui, J. Harting, and C. Misbah, "Two-dimensional vesicle dynamics under shear flow: Effect of confinement," *Phys. Rev. E* **83**, 066319 (2011).

19. D. Salac and M. Miksis, "Reynolds number effects on lipid vesicles," *J. Fluid Mech.* **711**, 122–146 (2012).

20. A. Laadhari, P. Saramito, C. Misbah, and G. Szekely, "Fully implicit methodology for the dynamics of biomembranes and capillary interfaces by combining the level set and newton methods," *J. Comput. Phys.* **343**, 271–299 (2017).

21. N. Valizadeh and T. Rabczuk, "Isogeometric analysis of hydrodynamics of vesicles using a monolithic phase-field approach," *Comput. Methods Appl. Mech. Eng.* **388**, 114191 (2022).

22. A. Laadhari, "Implicit finite element methodology for the numerical modeling of incompressible two-fluid flows with moving hyperelastic interface," *Appl. Math. Comput.* **333**, 376–400 (2018).

23. A. Laadhari, "An operator splitting strategy for fluid-structure interaction problems with thin elastic structures in an incompressible newtonian flow," *Appl. Math. Lett.* **81**, 35–43 (2018).

24. A. Torres-Sánchez, D. Millan, and M. Arroyo, "Modelling fluid deformable surfaces with an emphasis on biological interfaces," *J. Fluid Mech.* **872**, 218–271 (2019).

25. G. R. Cokelet, E. Merrill, E. Gilliland, H. Shin, A. Britten, and R. W. Jr, "The rheology of human blood—measurement near and at zero shear rate," *Transactions of the Society of Rheology* **7**, 303–317 (1963).

26. P. Neofytou and D. Drikakis, "Non-newtonian flow instability in a channel with a sudden expansion," *J. Non-Newton. Fluid Mech.* **111**, 127–150 (2003).

27. A. Copley, "Apparent viscosity and wall adherence of blood systems," (1960).

28. T. Willmore, *Riemannian geometry* (Oxford University Press, USA, 1993) ISBN = 0198514921.

29. F. Feng and W. S. Klug, "Finite element modeling of lipid bilayer membranes," *J. Comput. Phys.* **220**, 394–408 (2006).

30. A. Laadhari and G. Szekely, "Fully implicit finite element method for the modeling of free surface flows with surface tension effect," *Int. J. Numer. Methods Eng.* **111**, 1047–1074 (2017).

31. F. J. Walburn and D. J. Schneck, "A constitutive equation for whole human blood," *Biorheology* **13**, 201–210 (1976).

32. J. Janela, A. Lefebvre, and B. Maury, "A penalty method for the simulation of fluid - rigid body interaction," *ESAIM: Proc.* **14**, 115–123 (2005).

33. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, C. R. A. Logg, J. Ring, M. E. Rognes, and G. N. Wells, "The FEniCS project version 1.5," *Archive of Numerical Software* **3** (2015), 10.11588/ans.2015.100.20553.

34. R. J. LeVeque, "High-resolution conservative algorithms for advection in incompressible flow," *SIAM J. Numer. Anal.* **33**, 627–665 (1996).

35. H. Zhao and E. S. G. Shaqfeh, "The dynamics of a vesicle in simple shear flow," *J. Fluid Mech.* **674**, 578–604 (2011).

36. M. Kraus, W. Wintz, U. Seifert, and R. Lipowsky, "Fluid vesicles in shear flow," *Phys. Rev. Lett.* **77**, 3685 (1996).

37. V. Kantsler and V. Steinberg, "Orientation and dynamics of a vesicle in tank-treading motion in shear," *Phys. Rev. Lett.* **95**, 258101 (2005).

38. T. M. Fischer, M. Stöhr-Liesen, and H. Schmid-Schönbein, "The red cell as a fluid droplet: Tank tread-like motion of the human erythrocyte membrane in shear flow," *Science* **202**, 894–896 (1978).

39. A. Gizzi, R. Ruiz-Baier, S. Rossi, A. Laadhari, C. Cherubini, and S. Filippi, "A three-dimensional continuum model of active contraction in single cardiomyocytes," *Modeling, Simulation and Applications* **14**, 157–176 (2015).