## Abstract

In a suspension of extended objects such as colloidal particles, capsules or vesicles, the contribution of particles to the stress is usually evaluated by first determining the stress originating from a single particle (e.g. via integrating the fluid stress over the surface of a particle) and then adding up the contributions of individual particles. While adequate for a computation of the average stress over the entire system, this approach fails to correctly reproduce the *local* stress. In this work, we propose and validate a variant of the method of planes which overcomes this problem. The method is particularly suited for many-body interactions arising from, for example, shear and bending rigidity of red blood cells.

## 1. Introduction

Understanding the rheology of dense suspensions and emulsions is a complex field in modern physics and engineering. For example, it is of practical importance to find reliable constitutive laws for the rheology of suspensions of deformable particles, e.g. capsules, vesicles and red blood cells (RBCs), derived from microscopic properties such as particle shape and deformability. One quantity of central interest in such investigations is the stress tensor. The total stress is the sum of the fluid and particle stresses, ** σ**(

**,**

*r**t*)=

*σ*^{f}(

**,**

*r**t*)+

*σ*^{p}(

**,**

*r**t*). Usually, the evaluation of the local fluid stress

*σ*^{f}(

**,**

*r**t*) does not pose a serious problem [1–3]. The major difficulty is the computation of the local particle stress

*σ*^{p}(

**,**

*r**t*).

One may follow Batchelor [4] and use the fluid stress information at the surface of the particles or—equivalently—employ the approach presented by Ramanujan & Pozrikidis [5], who use the known membrane forces instead. When summed over all the particles, the correct volume-averaged particle stress may then be obtained. However, as will be shown in this paper, such ‘particle-based’ methods yield incorrect results on the *local* stress.

A common way to circumvent this problem is to study simple cases in which the macroscopic stress is constant over the entire channel cross section (e.g. planar Couette flow) or is known from analytical solutions of the momentum balance equation, as in the case of a Poiseuille flow (e.g. [6]) or Kolmogorov flow (e.g. [7]). A drawback of such a detour is that fluctuations of local shear stress are not accessible, thus setting severe limits to an analysis of the underlying problem within the framework of statistical physics (recall, for example, that equilibrium fluctuations of shear stress yield an independent measure of shear viscosity).

The present article focuses on this issue for the case of suspensions of soft objects confined between two planar walls. Here, the term ‘soft’ means that interaction potentials are differentiable so that a force can always be defined. We propose a method (a variant of the method of planes (MOP)) that allows the particle contribution to the local stress at any desired distance from a solid wall to be directly computed from interaction forces. As will be shown below, this method accurately reproduces the expected macroscopic behaviour of the local viscous shear stress while at the same time providing access to spatio-temporal fluctuations of this important quantity.

## 2. Stress evaluation

We study a dense suspension of RBCs immersed in a Newtonian fluid with shear viscosity *η*_{0}, which fills both the interior and exterior regions of the particles. Our simulation algorithm consists of a fluid solver, a membrane model and a fluid–membrane coupling scheme. For the fluid, we use the D3Q19 Bhatnagar–Gross–Krook lattice Boltzmann method (LBM) [8] with a linearized equilibrium function in order to avoid inertia [9]. The interaction of the fluid and the membrane is captured by the immersed boundary method (IBM) [10]. The membrane dynamics is derived from a constitutive model, and the strains in the membranes are evaluated using a finite-element method. A detailed description of the simulation method can be found in Krüger *et al.* [9].

As already mentioned, in a particulate suspension, the total stress is the sum of the fluid and particle stresses. Within the LBM, it is straightforward to evaluate the fluid stress locally (e.g. [2,3]). We thus focus here on a computation of the particle stress at a given position and independent of any further assumptions.

The effect of interactions on flow behaviour can be incorporated in the Navier–Stokes equations either (i) by direct implementation of particle forces as a (spatially and temporally varying) external force field or (ii) by introducing the particle stress tensor. The relation between membrane force density acting on the fluid and particle stress is
2.1In the following, we restrict the discussion to the case of a suspension confined between two parallel walls, separated by *L*_{z} along the *z*-direction. Noting that (*F*_{i} is the force on particle *i* and *r*_{i} is its position) and integrating equation (2.1) along the *x*- and *y*-directions, one arrives at
2.2where we defined , where *A*=*L*_{x}*L*_{y} (*L*_{x} and *L*_{y} being the lengths of the system along the *x*- and *y*-directions). Assuming that the particle contribution to the stress vanishes at the wall (which is always the case here), equation (2.2) can be integrated with respect to *z*. This gives
2.3where sgn is the sign function. In deriving equation (2.3), we used the fact that the volume integral over particle force is identically zero (momentum conservation). Equation (2.3) gives the stress on a plane at *z* parallel to the walls (figure 1).

The above approach is inspired by the so-called ‘MOP’, first introduced by Todd *et al.* [11] for the case of a simple liquid and further examined by Varnik *et al.* [12] in the case of a polymer melt. To the best of our knowledge, an extension of the MOP to the case of extended particles immersed in a fluid has not been proposed so far. The key point motivating the use of equation (2.3) in the case of interest to us is to realize that each individual membrane object (capsule, RBC) is represented by an ensemble of interacting nodes immersed in the ambient fluid. If we identify these membrane nodes as point particles, the fact that the thus defined point particles belong to extended objects plays no role in a computation of the particle stress. Thus, the index *i* in equation (2.3) specifies a node, and the sum runs over all the membrane nodes in the system. Note that the forces *F*_{i} in equation (2.3) do not necessarily have to be pair forces. This is essential for the applicability of the method in the present case since strain and bending resistances of the membranes effectively originate from many-particle interactions.

*A priori*, the fluid stress is known on the Eulerian lattice while forces needed for the evaluation of the particle stress are computed on the Lagrangian mesh. Owing to the fact that—in the present study—both systems are coupled by non-local IBM interpolations [9,10], this results in a weak but observable numerical noise of the total stress (data not shown). This problem is solved by first computing, via IBM interpolations, particle forces on the Eulerian grid and then applying the MOP to obtain the particle stress *Σ*^{p}_{xz}(*z*,*t*) completely in the Eulerian system. Equation (2.3) then reads
2.4where the sum runs over all lattice nodes *j* and *F*_{jx}(*t*) is the *x*-component of the force acting on node *j* at time *t*. This approach is reasonable since membrane forces are ‘visible’ through their action on the fluid only. As the fluid stress is computed at the fluid nodes, the particle stress should also be evaluated at those coordinates. For this purpose, we employ the relation , i.e. we take the average of two planes separated by one lattice constant *ΔZ*, where each plane is in the middle of two lattice grid lines (see figure 1). It is easy to see that this approach is equivalent to applying the MOP directly to a plane located at position *Z* on the Eulerian grid, ignoring all nodes on that plane. This is because, within the former approach, the contribution of the nodes directly on the plane at *Z* is taken twice but with different signs, thus leading to its cancellation.

## 3. Results

We consider the flow between two planar walls placed along the *z*-axis at positions *z*=0 and *z*=*D* and moving along the *x*-axis with velocities ±*u*_{w}, respectively (no external forces). This defines a global shear rate of . In such a set-up, it follows from the momentum balance equation that the *xz*-component (shear component) of the macroscopic deviatoric stress tensor is constant along the *z*-direction, *σ*_{xz}(*z*)=const. (‘macroscopic’ means that instantaneous fluctuations are averaged out). It is important to realize that this relation is valid for any flow profile and constitutive law of the fluid. It thus provides an excellent test case for our approach. In the following, all simulation quantities are given in lattice units, i.e. Δ*Z*=1, Δ*t*=1 and *ρ*=1.

We have performed six series of simulations with wall-to-wall separations *D*=15,…,90 in steps of 15. The size of the simulation boxes, including the bounce-back walls, is 90×90×(*D*+2) (table 1), with full periodicity along the *x*- and *y*-axes. The shear flow is imposed by moving bottom and top walls at *z*=0 and *z*=*D* via the bounce-back boundary condition [1]. In order to maintain the same average shear rate, the wall velocity increases linearly with the wall-to-wall distance and is 0.06 for the largest simulation box (table 1). Each series consists of 10 independent simulation runs with different random initial RBC configurations. The results in each series have been averaged in order to improve the ensemble averages. The LBM relaxation parameter is *τ*=1, and the two-point IBM interpolation stencil has been used [9]. The RBCs have a large radius *r*=9 and are discretized using 1620 triangular faces each (figure 2). The number of RBCs has been chosen in such a way that the haematocrit is 46.7 per cent in all studied cases (table 1). Each simulation ran for 50 000 time steps (according to nearly 67 inverse shear rates). The first 10 000 time steps have been omitted, thus neglecting transients. The remaining data have been used for time averaging (observables have been evaluated at every 100th time step).

Figure 3 contains the central result of our work. It compares the present method with a computation of local shear stress using ‘particle-based’ approaches, i.e. by computing the stress over each individual cell and then adding the contributions of all the cells whose centre of mass is within the interval [*Z*−1/2,*Z*+1/2]. While the latter approach fails to satisfy the requirement of a constant macroscopic shear stress across the channel, the proposed method perfectly fulfils this condition. We have seen this behaviour in all the cases studied. Generally, the ‘particle-based’ total stress reflects oscillations of the particle density (data not shown), while the new approach is free of this artefact. Indeed, such an ad hoc definition of the local stress must reflect, at least to some extent, the concentration profile of particles [12].

We emphasize that, when averaged over the entire simulation cell, both methods yield the same result—up to relative deviations of the order of 10^{−3} or less. The basic difference is thus in the evaluation of the local shear stress. This, however, is a major point when it comes to determining the local shear viscosity, which is the ratio of the local shear stress and the local rate of shear. As mentioned in §1, the knowledge of macroscopic stress may help to alleviate this problem. It does, however, not provide any access to the fluctuations of the local stress. In contrast to this, fluctuations of the local stress are a natural result of the present method (see the inset of figure 3*b*).

In figure 4, wall effects can easily be identified. In particular, the local viscosity is highest in a region close to the wall. The maximum and the location of this viscosity peak depend on the wall distance *D*. For a similar model to that in the present work, Doddi & Bagchi [6] also observed an increase in viscosity near the walls in a Poiseuille flow. Varnik & Binder [13], on the other hand, reported a significant increase (by two orders of magnitude) in the local viscosity of a polymer melt close to a solid substrate. However, suspensions of deformable particles usually exhibit a shear-thinning behaviour (e.g. the experiments of Chien [14] on RBCs), which depends on the applied shear rate or stress. In a Poiseuille flow, the stress varies linearly with distance from the wall and thus leads to a locally variable shear thinning [15], which superimposes the wall effects. Indeed, since the stress increases when approaching the walls, this would contribute to a corresponding *decrease* in viscosity owing to shear thinning. The observed increase in local shear viscosity in fig. 15 of Doddi & Bagchi [6] thus results from the competition between two opposite effects: the confinement effects—leading to an increase in local viscosity when approaching the wall—as opposed to variable stress effects that tend to decrease it. In the present study, on the other hand, the macroscopic stress is constant across the channel, which allows a separation of wall effects from influences related to shear thinning.

The strength of the present approach is further highlighted in figure 5. Clearly, the proposed method reproduces the expected constant macroscopic shear stress profile even in the case of a single capsule. Note that, in this case, the methods introduced by Batchelor [4], Ramanujan & Pozrikidis [5] and Ladd [1] would yield a single datum point as the average over the entire capsule.

## 4. Conclusions

A new approach for the evaluation of the local particle stress in immersed boundary lattice Boltzmann simulations is proposed and tested. It is possible to compute the particle stress as a function of the transverse coordinate independently of macroscopic assumptions. The spatial resolution is the same as that for the viscous fluid stress in the LBM.

The proposed method accurately recovers the constancy of shear stress in simple shear flow, even for a dense suspension of RBCs (46.7% volume fraction) with highly inhomogeneous viscosity. The volume average of the shear stress obtained with the newly proposed method equals the values calculated from two other methods which, however, fail to provide the correct behaviour of the local stress.

Even though there are ways, in the case of simple flow situations, to obtain the macroscopic average of the local shear stress, these approaches do not provide any access to stress fluctuations—a quantity of great importance within the framework of statistical physics. The problem is also solved within the present approach.

Finally, we show that the method works perfectly even in the case of a single object in a linear shear flow.

## Acknowledgements

This project has been supported by DFG grant VA205/5-1.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society