2.1. Previous models

Previous simulation studies, which have been performed to describe WBC margination, are based on strong model assumptions. The simulation approach (Freund Reference Freund2007) assumes an incompressible Stokes flow, the cells are modelled with a linear elastic membrane and a global area constraint is enforced. A boundary integral formulation is used for numerical discretization. More recently, a particle-based Lagrangian approach was used (Fedosov et al. Reference Fedosov, Fornleitner and Gompper2012; Fedosov & Gompper Reference Fedosov and Gompper2014). Here, RBCs and WBCs are described by a network model, where the cells are represented through triangulated surfaces. Penalty terms are used to ensure global volume and global area conservation as well as local area conservation for each surface element. The approach therefore guarantees inextensibility for sufficiently small surface elements. Each membrane point is connected to the fluid through viscous friction. The dynamics of the fluid flow is described by the smoothed dissipative particle dynamics (SDPD) method, an approximation for the Navier–Stokes equations which is only precise, if the particle density is large enough. Furthermore, the incompressibility of the fluid is not guaranteed a priori and has to be controlled. In Takeishi et al. (Reference Takeishi, Imai, Nakaaki, Yamaguchi and Ishikawa2014), a finite element approach is used for the RBCs, which are modelled as biconcave capsules and a Lattice–Boltzmann method for the fluid flow. The problems are coupled through an immersed boundary method.

Our approach is fully continuous and considers RBCs using a Helfrich-type curvature elastic model with an inextensibility constraint. WBCs are modelled either as rigid bodies or again with a Helfrich-type curvature elastic model, but with a constraint allowing for weak extensibility. Furthermore, we will consider the full Navier–Stokes equations to account for inertia effects in the plasma. According to previous studies (N’Dri, Shyy & Tran-Son-Tay Reference N’Dri, Shyy and Tran-Son-Tay2003) we also consider the internal fluid of the RBCs and WBCs as an incompressible Newtonian fluid, neglecting internal structures and the nucleus in WBCs.

2.2. Helfrich-type models

Helfrich-type modelling approaches have been applied to understand the complex motions and shape changes RBCs undergo within a flow field, e.g. tank-treating (TT) and tumbling (TB) motion (Fischer, Stöhr-Liesen & Schmid-Schönbein Reference Fischer, Stöhr-Liesen and Schmid-Schönbein1978). Within a low-Reynolds-number regime, the Stokes limit is valid and various numerical approaches have also been considered in this limit to analyse the TT and TB motion (Kraus et al. Reference Kraus, Wintz, Seifert and Lipowsky1996; Biben & Misbah Reference Biben and Misbah2003; Beaucourt et al. Reference Beaucourt, Rioual, Seon, Biben and Misbah2004; Biben, Kassner & Misbah Reference Biben, Kassner and Misbah2005; Veerapaneni et al. Reference Veerapaneni, Gueyffier, Zorin and Biros2009; Ghigliotti, Biben & Misbah Reference Ghigliotti, Biben and Misbah2010; Kim & Lai Reference Kim and Lai2010; Sohn et al. Reference Sohn, Tseng, Li, Voigt and Lowengrub2010; Zhao & Shaqfeh Reference Zhao and Shaqfeh2011). All models consider a membrane-free energy

(2.1) $$\begin{eqnarray}\displaystyle \mathscr{E}=\int _{{\it\Gamma}}\frac{1}{2}b_{N}(H-H_{0})^{2}\,\text{d}{\it\Gamma}, & & \displaystyle\end{eqnarray}$$

with membrane ${\it\Gamma}(t)$ , mean curvature $H$ , spontaneous curvature $H_{0}$ and normal bending rigidity $b_{N}$ . Lagrange multipliers are used to enforce a global area constraint or the stronger inextensibility constraint. The jump condition for the fluid stress tensor $\unicode[STIX]{x1D64E}=\unicode[STIX]{x1D61A}_{0,i}=-p\unicode[STIX]{x1D644}+{\it\nu}_{0,i}\unicode[STIX]{x1D63F}$ , with pressure $p$ , fluid viscosity ${\it\nu}_{0}$ , cell viscosity ${\it\nu}_{i}$ and deformation tensor $\unicode[STIX]{x1D63F}=\boldsymbol{{\rm\nabla}}\boldsymbol{v}+(\boldsymbol{{\rm\nabla}}\boldsymbol{v})^{\text{T}}$ , with velocity $\boldsymbol{v}$ , along the membrane then reads

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\boldsymbol{n}]_{{\it\Gamma}}=\frac{{\it\delta}\mathscr{E}}{{\it\delta}{\it\Gamma}}+{\it\lambda}_{\mathit{global}}\,H\boldsymbol{n}\quad \text{global area constraint,} & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\boldsymbol{n}]_{{\it\Gamma}}=\frac{{\it\delta}\mathscr{E}}{{\it\delta}{\it\Gamma}}+{\it\lambda}_{\mathit{local}}\,H\boldsymbol{n}+\boldsymbol{{\rm\nabla}}_{\!{\it\Gamma}}{\it\lambda}_{\mathit{local}}\quad \text{local inextensibility constraint}, & \displaystyle\end{eqnarray}$$

with outer normal $\boldsymbol{n}$ and the surface gradient $\boldsymbol{{\rm\nabla}}_{\!{\it\Gamma}}=\unicode[STIX]{x1D64B}\boldsymbol{{\rm\nabla}}$ , where $\unicode[STIX]{x1D64B}=\unicode[STIX]{x1D644}-\boldsymbol{n}\otimes \boldsymbol{n}$ denotes the projection operator. Index 0 thereby denotes the fluid phase and index $i$ the cell phase. The Lagrange multipliers ${\it\lambda}_{\mathit{global}}$ and ${\it\lambda}_{\mathit{local}}$ are functionals of the fluid velocity $\boldsymbol{v}$ and are obtained by requiring $\text{d}/\text{d}t\int _{{\it\Gamma}}\,\text{d}{\it\Gamma}=\int _{{\it\Gamma}}H\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}=0$ for the global area constraint and $\boldsymbol{{\rm\nabla}}_{\!{\it\Gamma}}\boldsymbol{\cdot }\boldsymbol{v}=0$ along ${\it\Gamma}(t)$ for the local inextensibility constraint. The jump condition for the velocity reads $[\boldsymbol{v}]_{{\it\Gamma}}=0$ .

The linearity of the Stokes problem allows for efficient decoupled algorithms to solve for the Lagrange multipliers (Biben & Misbah Reference Biben and Misbah2003; Sohn et al. Reference Sohn, Tseng, Li, Voigt and Lowengrub2010; Zhao & Shaqfeh Reference Zhao and Shaqfeh2011). However, in regimes with the Reynolds number being of order unity or higher, the Stokes limit is at least questionable. Modelling approaches, which consider inertia effects, have been introduced recently (Laadhari, Saramito & Misbah Reference Laadhari, Saramito and Misbah2012; Salac & Miksis Reference Salac and Miksis2012; Aland et al. Reference Aland, Egerer, Lowengrub and Voigt2014). All have found that the classical TB behaviour is no longer observed at moderate Reynolds numbers. As such, a suppression of TB motion could have far reaching consequences also for the interaction of RBCs and thus also the margination of WBCs. In the following, we will consider the full Navier–Stokes equations, which read inside and outside the RBCs

(2.4) $$\begin{eqnarray}\displaystyle & {\it\rho}(\partial _{t}\boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{v})-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}=0, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$

with density ${\it\rho}={\it\rho}_{0,i}$ . The global area constraint can be treated explicitly, which was e.g. used by Bonito, Nochetto & Pauletti (Reference Bonito, Nochetto and Pauletti2011) within a front tracking method, proposed by Du, Li & Liu (Reference Du, Li and Liu2007) and used by Haußer et al. (Reference Haußer, Li, Lowengrub, Marth, Rätz and Voigt2013) and Marth & Voigt (Reference Marth and Voigt2014) for a phase-field model and also considered in Salac & Miksis (Reference Salac and Miksis2011) for a level-set approach. The local inextensibility constraint is more delicate and leads to additional nonlinearities, which are so far only considered within a level-set approach in Salac & Miksis (Reference Salac and Miksis2011), Laadhari et al. (Reference Laadhari, Saramito and Misbah2012) and within a phase field approximation in Aland et al. (Reference Aland, Egerer, Lowengrub and Voigt2014). These models predominantly consider only one cell or simple test cases with more cells, using one level-set or phase field function for all cells. We will demonstrate, that this approach can lead to non-physical behaviour in more complex situations. We further like to distinguish between RBCs and WBCs having different properties, which can more naturally be considered by using different functions to describe different cells. It therefore will be the main modelling contribution to extent the described approaches in Aland et al. (Reference Aland, Egerer, Lowengrub and Voigt2014) to also model interactions of cells in an efficient way.

We will consider WBCs using two different approaches. The first models WBCs as a rigid body with a spherical shape using the fluid particle dynamics (FPD) approach (Tanaka & Araki Reference Tanaka and Araki2000), and the second accounts again for a Helfrich-type energy equation (2.1) with the jump condition for the fluid stress tensor described by

(2.6) $$\begin{eqnarray}\displaystyle [\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\boldsymbol{n}]_{{\it\Gamma}}=\frac{{\it\delta}\mathscr{E}}{{\it\delta}{\it\Gamma}}+\frac{c}{2}(\mathscr{A}^{0}-\mathscr{A}({\it\Gamma}))H\boldsymbol{n}\quad \text{weak global area constraint}, & & \displaystyle\end{eqnarray}$$

with penalty parameter $c$ and initial and desired cell area $\mathscr{A}^{0}$ and $\mathscr{A}({\it\Gamma})=\int _{{\it\Gamma}}\,\text{d}{\it\Gamma}$ , respectively. Appropriately choosing $c$ allows for moderate changes in area and thus an effective weak extensibility. The internal fluid of the WBCs are again considered as an incompressible Newtonian fluid, but with a larger viscosity as for the RBCs.

2.3. Hydrodynamic phase field models

The method introduces auxiliary phase fields ${\it\phi}_{i}$ that distinguish the inside and the outside of each cell $i=1,\ldots ,n$ . The inside and the outside are separated from each other by a diffuse layer, which marks the membrane. The phase field variables are defined as

(2.7) $$\begin{eqnarray}{\it\phi}_{i}(t,\boldsymbol{x}):=\tanh \left(\frac{r_{i}(t,\boldsymbol{x})}{\sqrt{2}{\it\epsilon}}\right)\!,\end{eqnarray}$$

where ${\it\epsilon}$ characterizes the thickness of the diffuse interface and $r_{i}(t,\boldsymbol{x})$ denotes the signed-distance function between $\boldsymbol{x}\in {\it\Omega}$ and its nearest point on ${\it\Gamma}_{i}(t)$ the membrane of cell $i$ . Depending on $r_{i}$ , we label the inside with ${\it\phi}_{i}\approx 1$ and the outside with ${\it\phi}_{i}\approx -1$ . ${\it\Gamma}_{i}(t)$ is then implicitly defined by the zero level set of ${\it\phi}_{i}$ . The cell phase can thus be defined as ${\it\phi}_{\mathit{cell}}\approx 1$ with ${\it\phi}_{\mathit{cell}}=\max _{\boldsymbol{x}\in {\it\Omega}}({\it\phi}_{1},\ldots ,{\it\phi}_{n})$ and the fluid phase as ${\it\phi}_{0}=-{\it\phi}_{\mathit{cell}}\approx 1$ . In the same way we define the RBCs phase ${\it\phi}_{RBCs}$ .

The dynamics are now governed by equations that couple these phase fields to the actual physical degrees of freedom. For WBCs, if considered as rigid bodies, the shape is a sphere and its velocity defined as the average fluid velocity inside

(2.8) $$\begin{eqnarray}\boldsymbol{v}_{i}(t)=\frac{\displaystyle \int _{{\it\Omega}}{\it\psi}_{i}(t,\boldsymbol{x})\boldsymbol{v}\,\text{d}{\it\Omega}}{\displaystyle \int _{{\it\Omega}}{\it\psi}_{i}(t,\boldsymbol{x})\,\text{d}{\it\Omega}},\end{eqnarray}$$

with ${\it\psi}_{i}(t,\boldsymbol{x})=0.5({\it\phi}_{i}(t,\boldsymbol{x})+1)$ . The motion of the WBC can then be described by $\boldsymbol{x}_{i}(t+{\rm\Delta}t)=\boldsymbol{x}_{i}(t)+{\rm\Delta}t\boldsymbol{v}_{i}(t)$ , with $\boldsymbol{x}_{i}(t)$ the centre of mass of the WBC and ${\rm\Delta}t$ the simulation time step (see Tanaka & Araki Reference Tanaka and Araki2000). In all other cases we consider the diffuse non-dimensional Helfrich energies

(2.9) $$\begin{eqnarray}\displaystyle \mathscr{E}_{i}({\it\phi}_{i})=\frac{1}{2\mathit{Re}\,\mathit{Be}_{i}}\int _{{\it\Omega}}\frac{1}{{\it\epsilon}}\left({\it\epsilon}{\rm\Delta}{\it\phi}_{i}-\frac{1}{{\it\epsilon}}({\it\phi}_{i}^{2}-1)({\it\phi}_{i}+H_{0})\right)^{2}\,\text{d}{\it\Omega}, & & \displaystyle\end{eqnarray}$$

with Reynolds number $\mathit{Re}={\it\rho}UL/{\it\nu}_{0}$ and bending capillary numbers $\mathit{Be}_{i}=(4\sqrt{2}/3)({\it\nu}_{0}UL^{2}/b_{N,i})$ , where $U$ denotes a characteristic velocity, $L$ a characteristic length and $b_{N,i}$ the bending rigidity of cell $i$ . In Du et al. (Reference Du, Liu, Ryham and Wang2005) formal convergence for ${\it\epsilon}\rightarrow 0$ to the non-dimensional form of the sharp interface energy in (2.1) is shown.

Instead of a direct extension of the models in Aland et al. (Reference Aland, Egerer, Lowengrub and Voigt2014), which considers an $L^{2}$ -gradient flow for ${\it\phi}_{i}$ and enforces a volume and local or global area constraint by Lagrange multipliers, we here introduce a $H^{-1}$ gradient flow. This directly ensures volume conservation and only area constraint have to be considered. They are treated as a global area constraint using a penalty approach. Such an approach was already considered for one cell, but without flow interactions, in Campelo & Hernández-Machado (Reference Campelo and Hernández-Machado2007). It has the advantage to keep the equations local. We introduce the non-dimensional penalty energies

(2.10) $$\begin{eqnarray}\displaystyle \mathscr{E}_{i,\mathit{area}}({\it\phi}_{i})=\frac{c_{i}}{2\mathit{Re}\,\mathit{Be}_{i}}(\mathscr{A}_{i}^{0}-\mathscr{A}({\it\phi}_{i}))^{2}, & & \displaystyle\end{eqnarray}$$

with penalty parameters $c_{i}$ , with different values for RBCs and WBCs to ensure global area conservation and weak extensibility, respectively. The initial and desired area of cell $i$ are denoted by $\mathscr{A}_{i}^{0}$ and $\mathscr{A}({\it\phi}_{i})=\int _{{\it\Omega}}({\it\epsilon}/2)|\boldsymbol{{\rm\nabla}}{\it\phi}_{i}|^{2}+(1/4{\it\epsilon})({\it\phi}_{i}^{2}-1)^{2}\,\text{d}{\it\Omega}$ , respectively. The last term converges to $(2\sqrt{2}/3)\int _{{\it\Gamma}_{i}}\,\text{d}{\it\Gamma}$ if ${\it\epsilon}\rightarrow 0$ (see Du, Liu & Wang Reference Du, Liu and Wang2006). The penalty energy therefore corresponds to the penalty term considered in (2.6). In case of an enforced stronger local inextensibility constraint for the RBCs, the penalty term helps to control the accumulation of errors, as shown in Du et al. (Reference Du, Liu and Wang2006) and Aland et al. (Reference Aland, Egerer, Lowengrub and Voigt2014).

In addition, we require an interaction energy $\mathscr{E}_{int}$ , to be defined below, such that the overall energy can be written as

(2.11) $$\begin{eqnarray}\mathscr{E}({\it\phi}_{1},\ldots ,{\it\phi}_{n})=\mathop{\sum }_{i=1}^{n}(\mathscr{E}_{i}({\it\phi}_{i})+\mathscr{E}_{i,area}({\it\phi}_{i}))+\mathscr{E}_{int}({\it\phi}_{1},\ldots ,{\it\phi}_{n})\end{eqnarray}$$

and the evolution equations for ${\it\phi}_{i}$ read

(2.12) $$\begin{eqnarray}\displaystyle \partial _{t}{\it\phi}_{i}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}_{i}={\it\gamma}{\rm\Delta}{\it\phi}_{i}^{\natural }, & & \displaystyle\end{eqnarray}$$

with a small positive mobility coefficient ${\it\gamma}$ and the non-dimensional chemical potentials

(2.13) $$\begin{eqnarray}\displaystyle {\it\phi}_{i}^{\natural }=\frac{{\it\delta}\mathscr{E}({\it\phi}_{1},\ldots ,{\it\phi}_{n})}{{\it\delta}{\it\phi}_{i}}=\frac{{\it\delta}\mathscr{E}_{i}({\it\phi}_{i})}{{\it\delta}{\it\phi}_{i}}+\frac{{\it\delta}\mathscr{E}_{i,area}({\it\phi}_{i})}{{\it\delta}{\it\phi}_{i}}+\frac{{\it\delta}\mathscr{E}_{int}({\it\phi}_{1},\ldots ,{\it\phi}_{n})}{{\it\delta}{\it\phi}_{i}}. & & \displaystyle\end{eqnarray}$$

We obtain

(2.14) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{{\it\delta}\mathscr{E}_{i}}{{\it\delta}{\it\phi}_{i}}=\frac{1}{\mathit{Re}\,\mathit{Be}_{i}}{\it\psi}_{i},\quad {\it\psi}_{i}={\rm\Delta}{\it\mu}_{i}-\frac{1}{{\it\epsilon}^{2}}(3{\it\phi}_{i}^{2}+2H_{0}{\it\phi}-1){\it\mu}_{i},\\ \displaystyle {\it\mu}_{i}={\it\epsilon}{\rm\Delta}{\it\phi}_{i}-\frac{1}{{\it\epsilon}}({\it\phi}_{i}^{2}-1)({\it\phi}_{i}+H_{0}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The penalty terms read

(2.15) $$\begin{eqnarray}\displaystyle \frac{{\it\delta}\mathscr{E}_{i,area}}{{\it\delta}{\it\phi}_{i}}=\frac{c_{i}}{\mathit{Re}\,\mathit{Be}_{i}}{\it\kappa}_{i}(\mathscr{A}_{i}^{0}-\mathscr{A}({\it\phi}_{i})),\quad {\it\kappa}_{i}={\it\epsilon}{\rm\Delta}{\it\phi}_{i}-\frac{1}{{\it\epsilon}}({\it\phi}_{i}^{2}-1){\it\phi}_{i}. & & \displaystyle\end{eqnarray}$$

Now, we have to consider the interaction terms. Interaction in principle is computationally costly, as it turns the problem into a nonlocal one and requires the coupling of all phase field variables ${\it\phi}_{1},\ldots ,{\it\phi}_{n}$ and computations of the distance between cells. We here consider only steric interactions to prevent coalescence or overlapping of cells and model the short range repulsion by a Gaussian potential, which in the sharp interface description reads

(2.16) $$\begin{eqnarray}\mathscr{E}_{i,int}({\it\Gamma}_{1},\ldots ,{\it\Gamma}_{n})=\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{n}{\it\alpha}\int _{{\it\Gamma}_{i}}w_{j}\,\text{d}{\it\Gamma},\quad \text{with }w_{j}(\boldsymbol{x})=\exp \left(-\frac{r_{j}^{2}(\boldsymbol{x})}{{\it\epsilon}^{2}}\right),\end{eqnarray}$$

where $w_{j}$ is an interaction function and describes the influence of cell $j$ on its environment. The interaction parameter ${\it\alpha}>0$ determines the strength of the repulsive interaction between cell $i$ and cell $j$ with respect to the evolution of cell $i$ . Using (2.7), the signed distance function $r_{j}$ can be computed within the diffuse interface region as

(2.17) $$\begin{eqnarray}r_{j}=-\frac{{\it\epsilon}}{\sqrt{2}}\ln \frac{1+{\it\phi}_{j}}{1-{\it\phi}_{j}}\quad \forall \boldsymbol{x}\boldsymbol{ : }|{\it\phi}_{j}(\boldsymbol{x})|<1.\end{eqnarray}$$

We thus can write the short-range interaction function $w_{j}$ as

(2.18) $$\begin{eqnarray}w_{j}=\left\{\begin{array}{@{}ll@{}}\exp \left(\displaystyle -\frac{1}{2}\left(\displaystyle \ln \frac{1+{\it\phi}_{j}}{1-{\it\phi}_{j}}\right)^{2}\right)\quad & \text{if }|{\it\phi}_{j}(\boldsymbol{x})|<1\\ 0\quad & \text{otherwise}\end{array}\right.\end{eqnarray}$$

and consider the interaction potential within the phase-field description, which reads in non-dimensional form

(2.19) $$\begin{eqnarray}\mathscr{E}_{i,int}({\it\phi}_{1},\ldots ,{\it\phi}_{n})=\frac{1}{\mathit{Re}\,\mathit{In}}\int _{{\it\Omega}}B({\it\phi}_{i})\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{n}w_{j}\,\text{d}{\it\Omega},\end{eqnarray}$$

with $B({\it\phi}_{i})=(1/{\it\epsilon})({\it\phi}_{i}^{2}-1)^{2}$ being non-zero only within the diffuse interface around ${\it\Gamma}_{i}$ and $\mathit{In}=(4\sqrt{2}/3)({\it\nu}_{0}U/{\it\alpha})$ the interaction number. The interaction energy thus reads

(2.20) $$\begin{eqnarray}\mathscr{E}_{int}({\it\phi}_{1},\ldots ,{\it\phi}_{n})=\mathop{\sum }_{i=1}^{n}\mathscr{E}_{i,int}({\it\phi}_{1},\ldots ,{\it\phi}_{n})\end{eqnarray}$$

and we obtain

(2.21) $$\begin{eqnarray}\displaystyle \frac{{\it\delta}\mathscr{E}_{int}({\it\phi}_{1},\ldots ,{\it\phi}_{n})}{{\it\delta}{\it\phi}_{i}}=\frac{1}{\mathit{Re}\,\mathit{In}}\left(B^{\prime }({\it\phi}_{i})\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{n}w_{j}+w_{i}^{\prime }\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{n}B({\it\phi}_{j})\right), & & \displaystyle\end{eqnarray}$$

with

(2.22) $$\begin{eqnarray}w_{i}^{\prime }=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{2}{{\it\phi}_{i}^{2}-1}\ln \frac{1+{\it\phi}_{i}}{1-{\it\phi}_{i}}\exp \left(-\frac{1}{2}\left(\ln \frac{1+{\it\phi}_{i}}{1-{\it\phi}_{i}}\right)^{2}\right)\quad & \text{if }|{\it\phi}_{i}(\boldsymbol{x})|<1\\ 0\quad & \text{otherwise}.\end{array}\right.\end{eqnarray}$$

Figure 1 gives a schematic illustration of the interaction terms. The algorithm considers only these cells, for which the diffuse interfaces overlap. All other cells do not contribute to the interaction. In addition, the most expensive part, computing the distance between cells, has been avoided, as this information is already contained in the phase field description of the cells. The approach thus scales with $n$ , the number of cells. Similar ideas to model interactions within phase field approaches have been considered in Zhang, Das & Du (Reference Zhang, Das and Du2009) and Gu, Wang & Gunzburger (Reference Gu, Wang and Gunzburger2014). However, only for the interaction of one cell with a fixed substrate.

Figure 1.

The red and blue coloured parts of ${\it\Gamma}_{i}$ are in contact with the interfaces of cell $j$ and cell $k$ (dashed contours around those cells), where the signed functions $r_{j}$ and $r_{k}$ can be calculated and thus also the interaction functions $w_{j}$ and $w_{k}$ . They do not vanish in the overlapping regions.

The non-dimensional Navier–Stokes equation reads

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}(\partial _{t}\boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{v})+\boldsymbol{{\rm\nabla}}p-\frac{1}{\mathit{Re}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\nu}\unicode[STIX]{x1D63F})=\mathop{\sum }_{i=1}^{n}{\it\phi}_{i}^{\sharp }\boldsymbol{{\rm\nabla}}{\it\phi}_{i}, & \displaystyle\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$

with ${\it\rho}=1$ and ${\it\nu}=(1-{\it\phi}_{\mathit{cell}}/2)+\sum _{i=1}^{n}({\it\nu}_{i}/{\it\nu}_{0})(({\it\phi}_{i}+1)/2)$ . Different densities could be handled in a similar way but are omitted here for simplicity.

To enforce the local inextensibility constraint for the RBCs, we follow the approach of model B in Aland et al. (Reference Aland, Egerer, Lowengrub and Voigt2014). The nonlinear evolution equations for ${\it\phi}_{i}$ remain, only the Navier–Stokes equation has to be extended and now reads

(2.25) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}(\partial _{t}\boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{v})+\boldsymbol{{\rm\nabla}}p-\frac{1}{\mathit{Re}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\nu}\unicode[STIX]{x1D63F})=\mathop{\sum }_{i=1}^{n}{\it\phi}_{i}^{\sharp }\boldsymbol{{\rm\nabla}}{\it\phi}_{i}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\frac{|\boldsymbol{{\rm\nabla}}{\it\phi}_{RBCs}|}{2}\unicode[STIX]{x1D64B}{\it\lambda}_{\mathit{local}}\right)\!,\quad & \displaystyle\end{eqnarray}$$

(2.26) $$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$

with a Lagrange multiplier ${\it\lambda}_{\mathit{local}}$ for which we introduce the additional equation

(2.27) $$\begin{eqnarray}{\it\xi}{\it\epsilon}^{2}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\phi}_{RBCs}^{2}\boldsymbol{{\rm\nabla}}{\it\lambda}_{\mathit{local}})+\frac{|\boldsymbol{{\rm\nabla}}{\it\phi}_{RBCs}|}{2}\unicode[STIX]{x1D64B}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{v}=0,\end{eqnarray}$$

with ${\it\xi}>0$ a parameter independent of ${\it\epsilon}$ . For ${\it\epsilon}\rightarrow 0$ we obtain ${\rm\Delta}{\it\lambda}_{\mathit{local}}=0$ away from ${\it\Gamma}_{RBCs}$ and $\unicode[STIX]{x1D64B}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{v}=\boldsymbol{{\rm\nabla}}_{\!{\it\Gamma}}\boldsymbol{\cdot }\boldsymbol{v}=0$ near ${\it\Gamma}_{RBCs}$ , which was shown in Aland et al. (Reference Aland, Egerer, Lowengrub and Voigt2014) for $n=1$ .

2.4. Thermodynamic consistency

The proposed system of equations (2.12) and (2.13) for $i=1,\ldots ,n$ and (2.23) and (2.24) or (2.25)–(2.27) fulfill thermodynamic consistency. To show this, we consider the kinetic energy $\mathscr{E}_{kin}=\int \boldsymbol{v}^{2}$ with constant density ${\it\rho}=1$ and the cell energy $\mathscr{E}$ and show that the time derivative is less than or equal to zero:

(2.28) $$\begin{eqnarray}\displaystyle \dot{\mathscr{E}}_{tot}(\boldsymbol{v},{\it\phi}_{1},\ldots ,{\it\phi}_{n})=\dot{\mathscr{E}}_{kin}+\dot{\mathscr{E}}=\int \boldsymbol{v}\boldsymbol{v}_{t}+\mathop{\sum }_{i=1}^{n}{\it\phi}_{i}^{\natural }\partial _{t}{\it\phi}_{i}\,\text{d}\boldsymbol{x}, & & \displaystyle\end{eqnarray}$$

with

(2.29) $$\begin{eqnarray}\displaystyle & \partial _{t}{\it\phi}_{i}=-\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}_{i}+{\it\gamma}{\rm\Delta}{\it\phi}_{i}^{\natural }, & \displaystyle\end{eqnarray}$$

(2.30) $$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}\boldsymbol{v}=-(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{v}-\boldsymbol{{\rm\nabla}}p+\frac{1}{\mathit{Re}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\nu}\unicode[STIX]{x1D63F})+\mathop{\sum }_{i=1}^{n}{\it\phi}_{i}^{\sharp }\boldsymbol{{\rm\nabla}}{\it\phi}_{i}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\frac{|\boldsymbol{{\rm\nabla}}{\it\phi}_{RBCs}|}{2}\unicode[STIX]{x1D64B}{\it\lambda}_{\mathit{local}}\right)\!,\qquad & \displaystyle\end{eqnarray}$$

which yields

(2.31) $$\begin{eqnarray}\displaystyle & & \displaystyle \dot{\mathscr{E}}_{tot}(\boldsymbol{v},{\it\phi}_{1},\ldots ,{\it\phi}_{n})\nonumber\\ \displaystyle & & \displaystyle \quad =\int \boldsymbol{v}\boldsymbol{\cdot }\left(-(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{v}-\boldsymbol{{\rm\nabla}}p+\frac{1}{\mathit{Re}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\nu}\unicode[STIX]{x1D63F})+\mathop{\sum }_{i=1}^{n}{\it\phi}_{i}^{\sharp }\boldsymbol{{\rm\nabla}}{\it\phi}_{i}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\frac{|\boldsymbol{{\rm\nabla}}{\it\phi}_{RBCs}|}{2}\unicode[STIX]{x1D64B}{\it\lambda}_{\mathit{local}}\right)\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\mathop{\sum }_{i=1}^{n}{\it\phi}_{i}^{\natural }(-\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}_{i}+{\it\gamma}{\rm\Delta}{\it\phi}_{i}^{\natural })\,\text{d}\boldsymbol{x}\quad (\text{partial integration, use }\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0)\nonumber\\ \displaystyle & & \displaystyle \quad =\int -\frac{1}{\mathit{Re}}|\boldsymbol{{\rm\nabla}}\boldsymbol{v}|^{2}-{\it\gamma}\mathop{\sum }_{i=1}^{n}|\boldsymbol{{\rm\nabla}}{\it\phi}_{i}^{\natural }|^{2}-\left(\boldsymbol{{\rm\nabla}}\boldsymbol{v}\boldsymbol{ : }\frac{|\boldsymbol{{\rm\nabla}}{\it\phi}_{RBCs}|}{2}\unicode[STIX]{x1D64B}\right){\it\lambda}_{\mathit{local}}\,\text{d}\boldsymbol{x}\nonumber\\ \displaystyle & & \displaystyle \qquad (\text{use (2.27) and partial integration})\nonumber\\ \displaystyle & & \displaystyle \quad =\int -\frac{1}{\mathit{Re}}|\boldsymbol{{\rm\nabla}}\boldsymbol{v}|^{2}-{\it\gamma}\mathop{\sum }_{i=1}^{n}|\boldsymbol{{\rm\nabla}}{\it\phi}_{i}^{\natural }|^{2}-{\it\xi}{\it\epsilon}^{2}{\it\phi}_{RBCs}^{2}|\boldsymbol{{\rm\nabla}}{\it\lambda}_{\mathit{local}}|^{2}\,\text{d}\boldsymbol{x}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant 0,\end{eqnarray}$$

where we have used the identity $\boldsymbol{v}\times (\boldsymbol{{\rm\nabla}}\times \boldsymbol{v})=\boldsymbol{{\rm\nabla}}(|\boldsymbol{v}|^{2})-(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{v}$ from which follows that $\int \boldsymbol{v}\boldsymbol{\cdot }(-\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{v})=0$ .

2.5. Numerical approach

Time discretization. In order to discretize in time, we explore an operator splitting approach. In an iterative process, we first solve the flow problem and substitute its solution into the phase-field equations, which are then solved separately with a parallel splitting method. We split the time interval $I=[0,T]$ into equidistant time instants $0=t_{0}<t_{1}<\cdots \,$ and define the time steps ${\it\tau}:=t_{n+1}-t_{n}$ . Of course, adaptive time steps may also be used. We define the discrete time derivative $\text{d}_{t}\cdot ^{n+1}:=(\cdot ^{n+1}-\cdot ^{n})/{\it\tau}$ , where the upper index denotes the time step number and e.g.  $\boldsymbol{v}^{n}:=\boldsymbol{v}(t_{n})$ is the value of $\boldsymbol{v}$ at time $t_{n}$ . For each system, a semi-implicit time discretization is used, which together with an appropriate linearization of the involved nonlinear terms leads to a set of linear systems in each time step.

Space discretization. We apply the finite element method to discretize in space, where a $P^{2}/P^{1}$ Taylor–Hood element is used for the flow problem, all other quantities are discretized in space using $P^{2}$ elements. In each time step we solve the following.

1. (i) The flow problem for $\boldsymbol{v}^{n+1}$ , $p^{n+1}$ and ${\it\lambda}_{\mathit{local}}^{n+1}$

   (2.32) $$\begin{eqnarray}\displaystyle \text{d}_{t}\boldsymbol{v}^{n+1}+(\boldsymbol{v}^{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{v}^{n+1} & = & \displaystyle -\boldsymbol{{\rm\nabla}}p^{n+1}+\frac{1}{\mathit{Re}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\nu}^{n}\unicode[STIX]{x1D63F}^{n+1})\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{i=1}^{n}{{\it\phi}_{i}^{\natural }}^{n}\boldsymbol{{\rm\nabla}}{\it\phi}_{i}^{n}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\frac{|\boldsymbol{{\rm\nabla}}{\it\phi}_{RBCs}^{n}|}{2}\unicode[STIX]{x1D64B}^{n}{\it\lambda}_{\mathit{local}}^{n+1}\right)\!,\end{eqnarray}$$
   (2.33) $$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}^{n+1}=0, & \displaystyle\end{eqnarray}$$
   (2.34) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\xi}{\it\epsilon}^{2}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(({\it\phi}_{RBCs}^{n})^{2}\boldsymbol{{\rm\nabla}}{\it\lambda}_{\mathit{local}}^{n+1})+\frac{|\boldsymbol{{\rm\nabla}}{\it\phi}_{RBCs}^{n}|}{2}\unicode[STIX]{x1D64B}^{n}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{v}^{n+1}=0, & \displaystyle\end{eqnarray}$$
   where ${\it\nu}^{n}={\it\nu}({\it\phi}^{n})$ and $\unicode[STIX]{x1D64B}^{n}=\unicode[STIX]{x1D644}-(\boldsymbol{{\rm\nabla}}{\it\phi}^{n}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}^{n}/|\boldsymbol{{\rm\nabla}}{\it\phi}^{n}|^{2})$ .

2. (ii) The phase field equations for ${\it\phi}_{i}^{n+1}$ , $i=1,\ldots ,n$

   (2.35) $$\begin{eqnarray}\text{d}_{t}{\it\phi}_{i}^{n+1}+\boldsymbol{v}^{n+1}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}_{i}^{n+1}={\it\gamma}{\rm\Delta}{{\it\phi}_{i}^{\natural }}^{n+1},\end{eqnarray}$$
   (2.36) $$\begin{eqnarray}\displaystyle {{\it\phi}_{i}^{\natural }}^{n+1} & = & \displaystyle \frac{1}{\mathit{Re}\,\mathit{Be}_{i}}{\it\psi}_{i}^{n+1}+\frac{c_{i}}{\mathit{Re}\,\mathit{Be}_{i}}(\mathscr{A}_{0,i}-\mathscr{A}({\it\phi}_{i}^{n})){\it\kappa}_{i}^{n}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{\mathit{Re}\,\mathit{In}}\left(B^{\prime }({\it\phi}_{i}^{n})\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{n}w_{j}^{n}+{w_{i}^{\prime }}^{n}\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{n}B({\it\phi}_{j}^{n})\right)\!,\end{eqnarray}$$
   (2.37) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\psi}_{i}^{n+1}={\rm\Delta}{\it\mu}_{i}^{n+1}-\frac{1}{{\it\epsilon}^{2}}(3({\it\phi}_{i}^{n+1})^{2}+2H_{0}{\it\phi}_{i}^{n+1}-1){\it\mu}_{i}^{n+1}, & \displaystyle\end{eqnarray}$$
   (2.38) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\mu}_{i}^{n+1}={\it\epsilon}{\rm\Delta}{\it\phi}_{i}^{n+1}-\frac{1}{{\it\epsilon}}(({\it\phi}_{i}^{n+1})^{2}-1)({\it\phi}_{i}^{n+1}+H_{0}), & \displaystyle\end{eqnarray}$$
   with ${\it\kappa}_{i}^{n}=-{\it\epsilon}{\rm\Delta}{\it\phi}_{i}^{n}+(1/{\it\epsilon})(({\it\phi}_{i}^{n})^{2}-1){\it\phi}_{i}^{n}$ . We linearize the nonlinear terms by a Taylor expansion of order one, e.g.  $(({\it\phi}_{i}^{n+1})^{2}-1){\it\phi}_{i}^{n+1}=(({\it\phi}_{i}^{n})^{2}-1){\it\phi}_{i}^{n}+(3({\it\phi}_{i}^{n})^{2}-1)$ $({\it\phi}_{i}^{n+1}-{\it\phi}_{i}^{n})$ .

Implementation. The fully discretized system of partial differential equations is implemented using the adaptive finite element toolbox AMDiS (Vey & Voigt Reference Vey and Voigt2007; Witkowski et al. Reference Witkowski, Ling, Praetorius and Voigt2015). We use an adaptively refined triangular mesh $\mathscr{T}_{h}$ with a high resolution along the cell membranes to guarantee at least five grid points across the diffuse interface. We further conduct a shared memory OPENMP parallelization, to solve the phase field evolutions via a parallel splitting method. Each linear system of equations is solved using the direct unsymmetric multifrontal method UMFPACK. For more details we refer to Ling et al. (Reference Ling, Marth, Praetorius and Voigt2015).