Lubricated Immersed Boundary Method in Two Dimensions

TL;DR

This paper introduces a lubricated immersed boundary method for 2D fluid-structure interaction problems involving thin fluid layers, achieving higher accuracy and better capturing vesicle dynamics near walls.

Contribution

The paper develops a novel immersed boundary method incorporating lubrication theory, improving accuracy in simulating near-contact elastic structures in fluid flows.

Findings

Achieves 2nd-order convergence in simple flows.

Effectively simulates vesicle near walls in shear flow.

Captures vesicle motion through narrow channels on coarse grids.

Abstract

Many biological examples of fluid-structure interaction, including the transit of red blood cells through the narrow slits in the spleen and the intracellular trafficking of vesicles into dendritic spines, involve the near-contact of elastic structures separated by thin layers of fluid. Motivated by such problems, we introduce an immersed boundary method that uses elements of lubrication theory to resolve thin fluid layers between immersed boundaries. We demonstrate 2nd-order accurate convergence for simple two-dimensional flows with known exact solutions to showcase the increased accuracy of this method compared to the standard immersed boundary method. Motivated by the phenomenon of wall-induced migration, we apply the lubricated immersed boundary method to simulate an elastic vesicle near a wall in shear flow. We also simulate the dynamics of a vesicle traveling through a narrow channel and observe the ability of the lubricated method to capture the vesicle motion on relatively coarse fluid grids.

Tables

Table 1. (b) L 1 superscript L 1 \textrm{L}^{1} errors and relative error in γ ˙ ˙ 𝛾 \dot{\gamma} .

gridsize	$u$	$v$	$p$	$𝐗$	$𝐔$	$\dot{\gamma }$
$16$	$4.4\cdot {10}^{-2}$	$2.2\cdot {10}^{-13}$	$1.3\cdot {10}^{-13}$	$4.7\cdot {10}^{-3}$	$3.0\cdot {10}^{-5}$	$1.5\cdot {10}^{-3}$
$32$	$1.8\cdot {10}^{-2}$	$8.0\cdot {10}^{-14}$	$4.9\cdot {10}^{-14}$	$5.3\cdot {10}^{-4}$	$1.5\cdot {10}^{-6}$	$7.7\cdot {10}^{-5}$
$64$	$6.3\cdot {10}^{-3}$	$5.8\cdot {10}^{-14}$	$2.0\cdot {10}^{-13}$	$6.6\cdot {10}^{-5}$	$7.6\cdot {10}^{-8}$	$3.9\cdot {10}^{-6}$
$128$	$2.3\cdot {10}^{-3}$	$1.8\cdot {10}^{-13}$	$2.8\cdot {10}^{-13}$	$1.1\cdot {10}^{-5}$	$4.7\cdot {10}^{-9}$	$2.4\cdot {10}^{-7}$

Table 2. (b) L 1 superscript L 1 \textrm{L}^{1} errors and relative error in γ ˙ ˙ 𝛾 \dot{\gamma} .

gridsize	$u$	$v$	$p$	$𝐗$	$𝐔$	$\dot{\gamma }$
$16$	$8.5\cdot {10}^{-4}$	$1.4\cdot {10}^{-14}$	$1.4\cdot {10}^{-14}$	$4.8\cdot {10}^{-4}$	$1.9\cdot {10}^{-8}$	$9.6\cdot {10}^{-7}$
$32$	$5.0\cdot {10}^{-4}$	$1.0\cdot {10}^{-14}$	$4.5\cdot {10}^{-15}$	$1.2\cdot {10}^{-4}$	$3.6\cdot {10}^{-9}$	$1.8\cdot {10}^{-7}$
$64$	$1.3\cdot {10}^{-4}$	$1.9\cdot {10}^{-14}$	$2.0\cdot {10}^{-13}$	$3.1\cdot {10}^{-5}$	$1.8\cdot {10}^{-10}$	$9.4\cdot {10}^{-9}$
$128$	$1.4\cdot {10}^{-4}$	$1.0\cdot {10}^{-13}$	$1.9\cdot {10}^{-13}$	$7.5\cdot {10}^{-6}$	$1.8\cdot {10}^{-10}$	$8.9\cdot {10}^{-9}$

Table 3. (a) L ∞ superscript L \textrm{L}^{\infty} errors

gridsize	$u$	$v$	$dp/d\theta$	$𝐗$	$𝐔$
$16$	$5.1\cdot {10}^{-3}$	$8.2\cdot {10}^{-3}$	$3.6\cdot {10}^{-1}$	$2.0\cdot {10}^{-3}$	$1.2\cdot {10}^{-5}$
$32$	$2.4\cdot {10}^{-3}$	$5.2\cdot {10}^{-3}$	$3.5\cdot {10}^{-1}$	$4.1\cdot {10}^{-4}$	$5.5\cdot {10}^{-7}$
$64$	$1.2\cdot {10}^{-3}$	$3.0\cdot {10}^{-3}$	$3.2\cdot {10}^{-1}$	$1.2\cdot {10}^{-4}$	$3.9\cdot {10}^{-7}$
$128$	$5.5\cdot {10}^{-4}$	$1.2\cdot {10}^{-3}$	$3.2\cdot {10}^{-1}$	$3.0\cdot {10}^{-5}$	$8.5\cdot {10}^{-8}$

Table 4. (b) L 1 superscript L 1 \textrm{L}^{1} errors

gridsize	$u$	$v$	$dp/d\theta$	$𝐗$	$𝐔$	$\dot{\gamma }$
$16$	$6.5\cdot {10}^{-3}$	$7.3\cdot {10}^{-3}$	$6.5\cdot {10}^{-2}$	$1.0\cdot {10}^{-2}$	$1.6\cdot {10}^{-5}$	$4.7\cdot {10}^{-3}$
$32$	$2.8\cdot {10}^{-3}$	$3.8\cdot {10}^{-3}$	$6.1\cdot {10}^{-2}$	$2.4\cdot {10}^{-3}$	$1.9\cdot {10}^{-6}$	$2.4\cdot {10}^{-5}$
$64$	$1.3\cdot {10}^{-3}$	$1.9\cdot {10}^{-3}$	$5.3\cdot {10}^{-2}$	$6.0\cdot {10}^{-4}$	$5.6\cdot {10}^{-7}$	$9.7\cdot {10}^{-6}$
$128$	$5.2\cdot {10}^{-4}$	$8.3\cdot {10}^{-4}$	$5.5\cdot {10}^{-2}$	$1.5\cdot {10}^{-4}$	$1.3\cdot {10}^{-7}$	$1.3\cdot {10}^{-6}$

Table 5. (a) L ∞ superscript L \textrm{L}^{\infty} errors

gridsize	$u$	$v$	$dp/d\theta$	$𝐗$	$𝐔$
$16$	$2.4\cdot {10}^{-3}$	$3.0\cdot {10}^{-3}$	$3.1\cdot {10}^{-1}$	$1.7\cdot {10}^{-3}$	$1.4\cdot {10}^{-6}$
$32$	$4.8\cdot {10}^{-3}$	$4.4\cdot {10}^{-3}$	$2.3\cdot {10}^{-1}$	$5.8\cdot {10}^{-4}$	$4.9\cdot {10}^{-7}$
$64$	$3.3\cdot {10}^{-3}$	$2.8\cdot {10}^{-3}$	$3.2\cdot {10}^{-2}$	$1.8\cdot {10}^{-4}$	$1.7\cdot {10}^{-7}$
$128$	$2.2\cdot {10}^{-3}$	$2.0\cdot {10}^{-3}$	$2.2\cdot {10}^{-2}$	$4.5\cdot {10}^{-5}$	$4.3\cdot {10}^{-8}$

Table 6. (b) L 1 superscript L 1 \textrm{L}^{1} errors

gridsize	$u$	$v$	$dp/d\theta$	$𝐗$	$𝐔$	$\dot{\gamma }$
$16$	$1.8\cdot {10}^{-3}$	$1.9\cdot {10}^{-3}$	$5.0\cdot {10}^{-2}$	$9.1\cdot {10}^{-3}$	$7.7\cdot {10}^{-6}$	$5.4\cdot {10}^{-6}$
$32$	$1.3\cdot {10}^{-3}$	$1.5\cdot {10}^{-3}$	$3.2\cdot {10}^{-2}$	$2.4\cdot {10}^{-3}$	$2.0\cdot {10}^{-6}$	$4.9\cdot {10}^{-7}$
$64$	$4.2\cdot {10}^{-4}$	$4.6\cdot {10}^{-4}$	$4.1\cdot {10}^{-3}$	$5.9\cdot {10}^{-4}$	$5.0\cdot {10}^{-7}$	$1.3\cdot {10}^{-7}$
$128$	$2.1\cdot {10}^{-4}$	$2.4\cdot {10}^{-4}$	$3.8\cdot {10}^{-3}$	$1.4\cdot {10}^{-4}$	$1.2\cdot {10}^{-7}$	$5.9\cdot {10}^{-8}$

Equations

$$\rho \frac{\partial \mathbf{u}}{\partial t}+\nabla p$$

$$\nabla \cdot \mathbf{u}$$

$$\mathbf{f}(\mathbf{x},t)={\int }_{\mathbf{q}}\mathbf{F}(\mathbf{q},t)\delta (\mathbf{x}-\mathbf{X}(\mathbf{q},t))\mathrm{d}\mathbf{q}.$$

$${\mathbf{U}}^{\text{IB}}(\mathbf{q},t)={\int }_{\Omega }\mathbf{u}(\mathbf{x},t)\delta (\mathbf{x}-\mathbf{X}(\mathbf{q},t))\mathrm{d}\mathbf{x}.$$

$${\mathbf{U}}_h+{\mathbf{U}}_0={\mathbf{U}}_h^{\text{IB}}+{\mathbf{U}}_0^{\text{IB}}$$

$$($$

$$\displaystyle\qquad\qquad\qquad\qquad\qquad+{\color[rgb]{0,0,0}\left((\mathbf{N_{0}}\cdot\nabla)\mathbf{u}\right)}\Big{|}_{0^{-}}\mkern-18.0mu\cdot\mathbf{T}_{0}+{\color[rgb]{0,0,0}\left((\mathbf{N_{0}}\cdot\nabla)\mathbf{u}\right)}\Big{|}_{h^{+}}\mkern-18.0mu\cdot\mathbf{T}_{0}\Bigg{)}$$

$$($$

$$\displaystyle\qquad\qquad\qquad-{\color[rgb]{0,0,0}(\mathbf{T_{0}}\cdot\nabla)}\left(h\mathbf{T}_{0}\cdot\left(\frac{\mathbf{U}_{h}+\mathbf{U}_{0}}{2}\right)\right)+\left(\mathbf{U}_{h}\cdot\mathbf{T}_{0}\right)({\color[rgb]{0,0,0}(\mathbf{T_{0}}\cdot\nabla)}h),$$

$$U_h-U_0={\int }_0^h\frac{\partial u}{\partial y}\mathrm{d}y,$$

$$U_{h}-U_{0}=h\frac{\partial u}{\partial y}\Bigg{|}_{y^{*}},$$

$$\frac{\partial p}{\partial x}$$

$$\frac{\partial p}{\partial y}$$

$$\partial u/\partial x+\partial v/\partial y$$

$$\frac{\partial u}{\partial y}=\frac{1}{2\mu }\frac{\partial p}{\partial x}(2y-h)+\frac{U_h-U_0}{h}.$$

$$\displaystyle U_{h}-U_{0}=\frac{h}{2}\left(\frac{\partial u}{\partial y}\Big{|}_{0^{+}}+\frac{\partial u}{\partial y}\Big{|}_{h^{-}}\right).$$

$$\displaystyle\frac{\partial p}{\partial x}=\frac{\mu}{h}\left(\frac{\partial u}{\partial y}\Big{|}_{h^{-}}-\frac{\partial u}{\partial y}\Big{|}_{0^{+}}\right).$$

$$u(y)=\frac{1}{2\mu }\frac{\partial p}{\partial x}y(y-h)+\frac{U_h-U_0}{h}y+U_0,$$

$$[{\sigma }_{ij}]N_j+F_i/J=0,$$

$${\sigma }_{ij}=-p{\delta }_{ij}+\mu \left(\partial u_i/\partial x_j+\partial u_j/\partial x_i\right).$$

$$-[p]N_i+\mu \left[\partial u_i/\partial x_j\right]N_j+\mu \left[\partial u_j/\partial x_i\right]N_j+F_j/J=0,$$

$$[p]\mathbf{N}-\mu \left[(\mathbf{N}\cdot \nabla )\mathbf{u}\right]=\mathbf{F}/J.$$

$$-\mu \left[(\mathbf{N}\cdot \nabla )\mathbf{u}\right]$$

$$[p]$$

$$\displaystyle\frac{\partial u}{\partial y}\Big{|}_{0^{+}}$$

$${\mathbf{N}}_0=({\mathbf{N}}_0\cdot {\mathbf{N}}_h){\mathbf{N}}_h+({\mathbf{N}}_0\cdot {\mathbf{T}}_h){\mathbf{T}}_h,$$

$$\left[\partial \mathbf{u}/\partial y\right]=\left[({\mathbf{N}}_0\cdot \nabla )\mathbf{u}\right]$$

$$=({\mathbf{N}}_0\cdot {\mathbf{N}}_h)\left[({\mathbf{N}}_{\mathbf{h}}\cdot \nabla )\mathbf{u}\right]$$

$$=({\mathbf{N}}_0\cdot {\mathbf{N}}_h)\left({\mathbf{F}}_h-({\mathbf{F}}_h\cdot {\mathbf{N}}_h){\mathbf{N}}_h\right)/J_h$$

$$\displaystyle\frac{\partial u}{\partial y}\Big{|}_{h^{-}}$$

$$U_{h}-U_{0}=\frac{h}{2}\left(\frac{(\mathbf{N}_{0}\cdot\mathbf{N}_{h})^{2}{\color[rgb]{0,0,0}(F_{h}^{\parallel}/J_{h})}-{\color[rgb]{0,0,0}(F_{0}^{\parallel}/J_{0})}}{\mu}+\frac{\partial u}{\partial y}\Bigg{|}_{0^{-}}+\frac{\partial u}{\partial y}\Bigg{|}_{h^{+}}\right),$$

$$V_h-V_0$$

$$=-{\int }_0^h\frac{\partial u}{\partial x}\mathrm{d}y$$

$$=-\left(\frac{\partial q}{\partial x}-U_h\frac{\partial h}{\partial x}\right),$$

$$q=\frac{-h^3}{12\mu }\frac{\partial p}{\partial x}+h\frac{U_h+U_0}{2},$$

$$V_h-V_0={\left(\frac{h^3}{12\mu }\frac{\partial p}{\partial x}\right)}_x-{\left(h\left(\frac{U_h+U_0}{2}\right)\right)}_x+U_h\frac{\partial h}{\partial x}.$$

Full text

††journal: Journal of Computational Physics

\stackMath

1 Introduction

The immersed boundary method is a widely-used numerical method for fluid-structure interaction that has been applied to problems including blood clotting [1], osmotic swelling due to thermal fluctuations [2], sperm locomotion through viscoelastic fluids [3], and insect flight [4]. Several extensions of the immersed boundary method have been developed to incorporate realistic structural properties such as added mass [5], permeability [6, 7], intrinsic twist [8], nonlinear constitutive laws [9], and variable viscosity and density [10, 11]. The chief idea of the immersed boundary method is to use an Eulerian description of the fluid and a Lagrangian description of the structure, while coupling these descriptions through integral operators involving delta functions. In the continuum formulation, these are the usual Dirac delta functions, and in practice much of the efficiency and accuracy of the immersed boundary method depends on how these Dirac delta functions are regularized [12, 13].

One of the features of the immersed boundary method is that it is formulated in a manner that avoids the contact problem; since all structures move in the same global velocity field, they cannot in principle cross themselves or one another [14]. Though it is useful for many purposes, advecting all structures in the same regularized velocity field presents challenges when multiple immersed boundaries are in near contact, which is typical in biological phenomena such as red blood cell motion through the microcirculation [15] and membrane receptor trafficking in dendrites [16]. In simulations using fixed fluid grids, the flow in the thin fluid layer between immersed boundaries may not be sufficiently resolved, resulting in boundaries that are effectively stuck together. This is sometimes dealt with by the addition of repulsive forces, as done by Krishnan et al. to prevent “unphysical overlap” of rigid particles [17], and by Lim et al. in their simulations of flagellar bundling in E. Coli to “prevent filaments from crossing each other” [18]. We wish to overcome this difficulty by directly using the governing equations instead.

The immersed boundary method takes the viewpoint that the entire domain, including any immersed elastic structures, is filled with fluid having a continuously varying velocity field, without large jumps in the velocities of structures that are nearby relative to the grid spacing $\Delta x$.

Given the difficulties observed above, we propose a modification to the immersed boundary method that makes use of lubrication theory to resolve the flow in thin layers. Unlike the standard immersed boundary method, which becomes less accurate as two boundaries approach each other, the lubrication approximation actually improves as the gap size decreases. The grid spacing $\Delta x$ provides a natural cutoff at which to apply the lubrication approximation. In the proposed method, the lubrication approximation is automatically turned off when the gap size becomes larger than the grid spacing, as described in detail later on.

We illustrate this method by applying it to two test problems of increasing complexity for which the lubricated immersed boundary method has significantly smaller errors than the standard method. The primary additional cost of the lubricated immersed boundary method involves the construction of a height function that gives the distance between the immersed boundaries. We describe a method to compute this height function accurately and efficiently. Finally, we apply the lubricated immersed boundary method to study the motion of elastic vesicles near walls and observe the significantly increased accuracy of simulations of wall-induced migration and channel flow on coarse fluid grids.

Our work builds on several previous efforts to merge direct numerical simulation with lubrication theory. In the context of particle suspensions, Nguyen and Ladd and Janoschek et al. worked out lubrication corrections to lattice-Boltzmann simulations of particle suspensions [19, 20], while Seto et al. performed simulations of hard spheres with regularized lubrication forces [21]. Thomas et al. performed a multiscale computation of a fluid drop interacting with a wall in which lubrication theory was used to replace the no-slip boundary condition by a specified wall-shear stress [22]. Our approach is most closely related to this last work, since we also combine a grid-based fluid solver with analytical results from lubrication theory to model subgrid effects. Moreover, there are common features in terms of the mathematical formulation; we also use a height function to describe the size of the lubricating layer between boundaries, and relate the flow in the gap to the stress on the boundaries. However, whereas Thomas et al. solve an evolution equation for the height and flux through the gap, the solution to which gives a boundary condition for the fluid, our boundaries are defined in a Lagrangian manner with the height function constructed based on an explicit representation of the boundaries as piecewise polynomial curves. In our approach, we use lubrication theory to modify the advection velocity. Further, whereas Thomas et al. considered a multiphase fluid in which a drop of fluid having no surface elasticity and with ten times greater density and viscosity than the ambient fluid falls down a slope, in the present study we restrict attention to elastic structures immersed in fluids with constant density and viscosity.

We note that an alternative approach for ensuring that the fluid dynamics is sufficiently resolved throughout the domain is to use adaptive mesh refinement (AMR). The idea of AMR is to use a hierarchy of fluid grids with different levels of refinement, with the finer grids placed in specific regions of space to obtain a desired accuracy [23]. This allows for highly-resolved computations at drastically reduced computational cost compared to uniform grids. However, implementing AMR requires additional data structures to keep track of the grid hierarchy, and typical implementations of AMR require the user to specify a finest grid level so that issues with near-contact are not definitively resolved. This finest grid size can also lead to severe timestep restrictions for numerical stability. The advantage of using a subgrid model, as we do here through lubrication theory, is that the asymptotic results hold for infinitesimally small gaps, and does not introduce additional timestep constraints. Further, relatively little additional machinery is required beyond what is already needed for uniform grid simulations; as we will discuss later on, the main extra cost of the lubricated immersed boundary method is the construction of the height function that gives the distance between the immersed boundaries.

The movies referenced throughout the manuscript are available at author TGF’s website at http://scholar.harvard.edu/tfai/home.

2 Formulation

First, we review the standard immersed boundary method, so that the new features of the lubricated method will be clear when they are introduced later on.

2.1 Classical Immersed Boundary Method

Consider a domain $\Omega$ filled with incompressible fluid of density $\rho$ and viscosity $\mu$, having velocity $\mathbf{u}(\mathbf{x},t)$ and pressure $p(\mathbf{x},t)$ defined in terms of cartesian coordinates $\mathbf{x}$. For many relevant problems, the nonlinear term in the Navier-Stokes equations may be neglected, and we will restrict attention to this case in which the fluid is described by the unsteady Stokes equations:

[TABLE]

where $\mathbf{f}(\mathbf{x},t)$ is the applied force density on the fluid.We are interested in the case in which an immersed elastic body exerts forces on the fluid and thereby determines $\mathbf{f}$.

The elastic body is described using Lagrangian coordinates $\mathbf{q}$, e.g. arclength in the undeformed frame for a one-dimensional immersed boundary, and the function $\mathbf{X}(\mathbf{q},t)$ gives the cartesian position at time $t$ of the material point labeled by $\mathbf{q}$ along the boundary. We assume that the Lagrangian force density $\mathbf{F}(\mathbf{q},t)$ is calculated in terms of $\mathbf{X}(\mathbf{q},t)$; for instance, if there exists an energy functional $E=E[\mathbf{X}(\cdot,t)]$, $\mathbf{F}$ may be calculated via $\mathbf{F}=-\varrho E/\varrho\mathbf{X}$, where $\varrho/\varrho\mathbf{X}$ denotes the first variation or Fréchet derivative. The Dirac delta function is used to convert the Eulerian force density from Lagrangian to cartesian coordinates via

[TABLE]

The immersed boundary moves at a velocity $\mathbf{U}^{\text{IB}}(\mathbf{q},t)$ equal to the local fluid velocity. This condition can be expressed in terms of the Dirac delta function as follows:

[TABLE]

As noted in [12], the above equation can be simplified using the definition of the Dirac delta function to obtain $\mathbf{U}^{\text{IB}}(\mathbf{q},t)=\mathbf{u}(\mathbf{X}(\mathbf{q},t),t)$, but the expression (4) is particularly convenient since this integral form is the basis of the numerical method used, in which the Dirac delta functions are replaced by regularized delta functions and $\mathbf{U}^{\text{IB}}(\mathbf{q},t)$ is evaluated by summation over the fluid grid. Regularization gives the immersed boundaries an effective width or blurriness on the order of the grid spacing $\Delta x$.

2.2 Lubrication Corrections

Now, consider two immersed boundaries separated by a thin layer of fluid with height $h\ll 1$. Note that throughout this article $h$ will be used to represent height and should not be assumed equal to the grid spacing $\Delta x$ (or $\Delta y$). In the limit $h\to 0$, the two immersed boundaries coincide and move at identical velocities given by (4). By the same logic, we expect that when boundary velocities are interpolated using regularized delta functions, boundaries separated by a gap smaller than the grid spacing $\Delta x$ will move at similar velocities and stick together. A thought experiment illustrates the issue: imagine two parallel lines separated by a distance $h$ and pulled in opposite directions. The Lagrangian forces on the lines are equal and opposite, and if the regularization length scale is large compared to the spacing between the lines, the delta function layers of force will be blurred so that they nearly cancel. In the limit $h\to 0$, this leads to a vanishing Eulerian force density in (3), and therefore vanishing velocity fields $\mathbf{u}(\mathbf{x},t)$ and $\mathbf{U}^{\text{IB}}(\mathbf{q},t)$. In practice, similar difficulties are encountered whenever $h\lesssim\Delta x$, since in this case the spread forces partially cancel and the thin lubricating layers separating immersed boundaries are not sufficiently resolved. Instead of computing the distinct velocities of the two boundaries, underresolved simulations assign the averaged velocity to both immersed boundaries, and the lines fail to move unless an extremely large force is exerted.

Motivated by this observation, we propose an alternative formulation for the advection velocity $\mathbf{U}(\mathbf{q},t)$ with the property $\mathbf{U}(\mathbf{q},t)\neq\mathbf{U}^{\text{IB}}(\mathbf{q},t)$. We use the average velocity computed by the standard immersed boundary method, but compute the velocity difference in an alternative way. Assume that there are two immersed boundaries close to touching, and let $\mathbf{U}_{0}$ and $\mathbf{U}_{h}$ be the velocities on the lower and upper surfaces, respectively, with corresponding force densities $\mathbf{F}_{0}$ and $\mathbf{F}_{h}$ and positions $\mathbf{X}_{0}$ and $\mathbf{X}_{h}$ satisfying $\mathbf{X}_{h}=\mathbf{X}_{0}+h\mathbf{N}_{0}$, where $\mathbf{N}_{0}$ is the unit normal to the lower surface at position $\mathbf{X}_{0}$ (Figure 1). In our formulation, $\mathbf{U}_{0}$ and $\mathbf{U}_{h}$ are determined by the following equations specifying their sum and difference:

[TABLE]

where $\mathbf{T}$ is the unit tangent vector and $J=\det\left(\partial X_{i}/\partial q_{j}\right)$ is the Jacobian, so that $\mathbf{F}/J$ has units of an Eulerian force density. We write $F^{\parallel}=\mathbf{T}\cdot\mathbf{F}$ and $F^{\perp}=\mathbf{N}\cdot\mathbf{F}$ to denote the tangential and normal components of the force density, respectively. The expressions $(\mathbf{N_{0}}\cdot\nabla)$ and $(\mathbf{T_{0}}\cdot\nabla)$ are used to denote directional derivatives normal and tangent to the lower surface, respectively. In the case of a flat lower boundary aligned with the $x$-axis, $(\mathbf{N_{0}}\cdot\nabla)\equiv\partial/\partial y$, and $(\mathbf{T_{0}}\cdot\nabla)\equiv\partial/\partial x$.

Together, (5) and (6)–(7) define $\mathbf{U}_{h}$ and $\mathbf{U}_{0}$; to summarize, the standard immersed boundary formulation is used for the average velocity and a sharp formulation is used for the velocity difference. It is reasonable to compute the mean velocity using the standard immersed boundary method as in (5), since smoothing tends to yield the net motion as in the thought experiment above. The distinction between our approach and the standard immersed boundary method is the formulation for the difference $\mathbf{U}_{h}-\mathbf{U}_{0}$, which has several practical advantages compared to the standard method. Consider (6), for instance: the Lagrangian forces $\mathbf{F}_{0}$ and $\mathbf{F}_{h}$ are quantities supported on the immersed boundary, and it follows that $F_{0}^{\parallel}$ and $F_{h}^{\parallel}$ are never smoothed by convolution against a delta function kernel. Moreover, the quantities $((\mathbf{N_{0}}\cdot\nabla)\mathbf{u})|_{0^{-}}$ and $\left((\mathbf{N_{0}}\cdot\nabla\right)\mathbf{u})|_{h^{+}}$ are outer derivatives, which implies they are typically well-resolved and/or small in magnitude. Therefore, the terms on the right-hand side of (6) are sharply defined and convenient to use, and (7) is similarly formulated in terms of quantities supported on the boundary or outside the lubrication layer. Next, we show how (6)–(7) may be derived from lubrication theory.

2.3 Lubrication Equations

In the thought experiment presented above involving two parallel lines pulled in opposite directions with $x$-velocities $U_{0}$ and $U_{h}$, the exact solution $\mathbf{u}=(u,v)$ inside the gap between the lines is linear shear flow, with $u(y)=(U_{h}-U_{0})y/h+U_{0}$ and $v\equiv 0$. We wish to use this formula for $u(y)$ to obtain the difference $U_{h}-U_{0}$ in terms of known quantities, and it will turn out that the derivative $\partial u/\partial y$ normal to the lines is a convenient choice. By the fundamental theorem of calculus,

[TABLE]

and because the derivative $\partial u/\partial y$ is constant for linear shear, it follows that

[TABLE]

where $y^{*}$ is any value of $y$ in $(0,h)$. Although the above equation is trivial, as it simply states that the flow in a linear shear is completely determined by its shear rate, surprisingly there is a similar statement for flows through narrow lubrication layers in general.

To prove this more general statement, we first review some classical results from lubrication theory, although we will use these results in a nonstandard manner. In most applications of lubrication theory, velocity boundary conditions are prescribed and the pressure is calculated by solving the Reynolds equation [24], whereas in our case the relative velocity of the boundaries is the unknown.

Assume a flat lower surface aligned with the $x$-direction is separated by height $h$ from a (possibly curved) upper surface with characteristic horizontal length scale $L$, where $h\ll L$. The unsteady Stokes equations can be simplified substantially in this case by neglecting terms of order $h/L$ and smaller to obtain [24]

[TABLE]

Since $p$ is $y$-independent, (10) may be integrated explicitly to yield

[TABLE]

Because the first term of the right-hand side has opposite signs for $y=0$ and $y=h$, we may eliminate the pressure gradient and invert the above equation for $U_{h}-U_{0}$ to get

[TABLE]

This implies that the observation (9) made in the special case of linear shear that the velocity difference may be computed by knowledge of $\partial u/\partial y$ at the boundaries applies to lubrication flows in general. In fact, the pressure through the gap can also be recovered using

[TABLE]

Knowledge of the pressure gradient and boundary velocities is sufficient to reconstruct the full tangential velocity profile through the gap by the formula

[TABLE]

which follows from integrating (13). (The normal velocities may be computed in terms of the pressure gradient and tangential velocities, as shown later on.)

We have now derived how the boundary velocities and flow through the gap may be calculated in terms of normal derivatives on the boundaries. Next, we will show these normal derivatives may themselves be calculated in terms of more convenient quantities.

2.4 Jump Conditions

Here, we review the classical jump conditions across an immersed boundary [25, 26], closely following the exposition of [27]. Let $[\phi]$ denote the jump in a quantity $\phi$ across an immersed boundary, so that e.g. $[\sigma_{ij}]$ for $i,\,j=1,\,2$ gives the jump in the components of the stress tensor. Since immersed boundaries are considered massless, all stresses must balance, i.e.

[TABLE]

where $\mathbf{F}$ is the Lagrangian force density, $\mathbf{N}$ is the unit normal vector, and $J$ is Jacobian. We will continue to assume for now that the lower surface is aligned with the $x$-axis so that $\mathbf{N}_{0}\equiv\mathbf{\hat{\j}}$. However, the upper surface is not assumed flat, i.e. $\mathbf{N}_{h}$ may vary from point to point depending on its curvature. As we will need to use the jump conditions across both surfaces, we use a more general notation that involves surface normals and tangents. In a Newtonian fluid the stress tensor is given by

[TABLE]

It follows that

[TABLE]

which can be simplified [27] to

[TABLE]

Separating (20) into tangential and normal components, we have

[TABLE]

where we have used that the normal component of $\left[(\mathbf{N}\cdot\nabla)\mathbf{u}\right]$ vanishes [27]. The first jump condition (21) relates forces on the immersed boundary to the normal derivative in velocity, which by (14) is just the quantity needed to find $U_{h}-U_{0}$. On the lower surface, $\mathbf{N}_{0}\equiv\mathbf{\hat{\j}}$ implies that $(\mathbf{N}_{0}\cdot\nabla)\mathbf{u}=\partial\mathbf{u}/\partial y$, and we may use (21) to rewrite the first term of (14) as

[TABLE]

where $F^{\parallel}=(\mathbf{F}\cdot\mathbf{T})$ is the tangential component of the Lagrangian force density and $\mathbf{T}$ is a unit tangent vector, i.e. $\mathbf{T}_{0}\equiv\mathbf{\hat{\i}}$. On the other hand, the unit normal $\mathbf{N}_{h}$ to the upper surface does not in general point in the $y$-direction, and so the associated jump condition does not directly involve $\partial\mathbf{u}/\partial y$. However, since

[TABLE]

we have for the jump across the upper surface that

[TABLE]

where we have used that there is no jump in tangential derivatives [27] in the second line, and the jump condition (21) in the final line. The right hand side above is clearly tangential to the upper surface, and therefore we may take the dot product with $\mathbf{T}_{0}$ to obtain $[\partial u/\partial y]$. Since $(\mathbf{T}_{0}\cdot\mathbf{T}_{h})=(\mathbf{N}_{0}\cdot\mathbf{N}_{h})$ in 2D this results in

[TABLE]

Substituting into (14) the expressions (23) and (26) yields

[TABLE]

To this point, we have not considered the normal velocities $V_{0}$ and $V_{h}$ on the lower and upper boundaries, respectively. In order to obtain an expression for $V_{h}-V_{0}$, we invoke the incompressibility constraint $\partial u/\partial x+\partial v/\partial y=0$ to write

[TABLE]

where $q=\int_{0}^{h}u\mathrm{d}y$ is the flow rate and we have used the fundamental theorem of calculus in the first and third equalities. Integrating (16) directly yields

[TABLE]

and substituting this expression into (28) results in

[TABLE]

Note the asymmetry with respect to the lower and upper surfaces, in that the final term on the right-hand side of (30) has no corresponding term involving $U_{0}$. This is a consequence of defining the height in terms of the lower surface normal.

As was the case for the tangential velocities, we would like to be able to compute $V_{h}-V_{0}$ in terms of convenient quantities such as the boundary forces. This is possible using the above formula (30), which involves tangential velocities already known through (27), the pressure gradient $\partial p/\partial x$ already known through (15), and derivatives of these quantities. In practice, to compute the pressure gradient we have found it preferable to use the jump condition (22), according to which

[TABLE]

where $F^{\perp}=(\mathbf{F}\cdot\mathbf{N})$. We average $p|_{0^{+}}$ and $p|_{h^{-}}$ to get a pressure in the gap and obtain the following formula:

[TABLE]

Thus we have an equation for $V_{h}-V_{0}$ that involves boundary forces, horizontal velocities that are themselves computed in terms of boundary forces, and the outer pressures $p|_{0^{-}}$ and $p|_{h^{+}}$.

For the general case of a nonlinear lower surface, we replace $x$- and $y$-velocities by velocities tangential and normal to the lower surface, respectively, and $x$-derivatives by tangential derivatives to the lower surface, resulting in the equations (6) and (7) for the lubricated immersed boundary method. Substituting coordinates in this manner is reasonable provided that $h\ll L$ so that the effects of curvature are negligible. The equations analogous to (27) and (33) in the general case are (6) and (7) of Section 2.2, respectively. The equation analogous to (16) for the tangential velocity through the gap in the case of curved surfaces is

[TABLE]

where $w\in[0,h]$ and where $\mathbf{x}=\mathbf{X}_{0}+w\mathbf{N}_{0}$ for fixed $\mathbf{X}_{0}$, so that $\mathbf{x}$ travels through the gap in the direction normal to the lower curve. Given this explicit representation for the tangential velocity, the normal component of the gap velocity may be computed as before via incompressibility:

[TABLE]

where the integrand can be computed through (34). Note that for the upper limit of integration in (35), $\mathbf{x}\cdot\mathbf{N}_{0}$ is simply the height above the lower surface of an arbitrary point $\mathbf{x}$ inside the gap.

3 Numerical Method

3.1 Spatial discretization

We next specialize to the case of a rectangular unit domain $\Omega$ of size $L_{x}$ by $L_{y}$ and periodic boundary conditions in both directions. We discretize $\Omega$ by overlaying an $N_{x}\times N_{y}$ grid with grid spacing $\Delta x=L_{x}/N_{x}$ and $\Delta=L_{y}/N_{y}$ in the $x$- and $y$-directions, respectively. The $(i,j)$th grid cell has its center at $\mathbf{x}_{i,j}=((i+1/2)\Delta x,(j+1/2)\Delta y)$ for $i=0,\dots,N_{x}-1$ and $j=0,\dots,N_{y}-1$.

We will consider a one-dimensional boundary immersed in $\Omega$ and parameterized by the Lagrangian coordinate $q\in[0,L_{q}]$. The immersed boundary is discretized using $N_{q}$ points so that the Lagrangian grid spacing is $\Delta q=L_{q}/N_{q}$, and the positions of Lagrangian points in cartesian coordinates are given by $\mathbf{X}_{k}(t):=\mathbf{X}(k\Delta q,t)$ for $k=0,\dots,N_{q}-1$.

The grid velocities $\mathbf{u}_{i,j}$ are located at cell centers $\mathbf{x}_{i,j}$, as are the external forces. The pressures are located at cell corners, with $p_{i,j}$ defined at the bottom-right corner of the $(i,j)$th grid cell. Our discretization is based on [28], which uses a finite element projection to obtain a symmetric positive definite system that can be solved efficiently by methods such as geometric multigrid. The unsteady Stokes equations are discretized in space in a centered, 2nd-order accurate manner:

[TABLE]

where

[TABLE]

and where all index arithmetic is to be interpreted modulo $N_{x}$ or $N_{y}$. We also have the equation of motion for the immersed boundary points

[TABLE]

To close the system, we must provide discretized equations for the spreading and interpolation operators used to compute $\mathbf{f}_{i,j}$ and $\mathbf{U}_{k}$. They are as follows:

[TABLE]

where we use the standard 4-point immersed boundary delta function $\delta_{h}(\mathbf{x})=\phi(x/\Delta x)\phi(y/\Delta y)$ [12], in which $\phi$ is given by

[TABLE]

We assume that the Lagrangian force $\mathbf{F}_{k}(t)$ can be calculated in terms of the positions $\mathbf{X}_{k}(t)$, specific examples of which will be given later on.

3.2 Timestepping

Upon discretizing in space, we have a semi-discrete problem with unknowns $\mathbf{u}_{i,j}(t)$, $p_{i,j}(t)$, and $\mathbf{X}_{k}(t)$ that depend on a continuous time variable. The problem is fully discretized by introducing a timestep $\Delta t$ and letting $\mathbf{u}^{n}_{i,j}:=\mathbf{u}_{i,j}(n\Delta t)$, $p^{n}_{i,j}:=p_{i,j}(n\Delta t)$, and $\mathbf{X}^{n}_{k}:=\mathbf{X}_{k}(n\Delta t)$. We use a second-order accurate predictor-corrector timestepping scheme. A preliminary approximation $\mathbf{X}^{*}$ to $\mathbf{X}^{n+1}$ is first computed through a first-order accurate forward Euler step:

[TABLE]

where $\mathbf{U}^{n}$ is computed in terms of the grid velocity $\mathbf{u}^{n}$ via (42). We calculate $\mathbf{u}^{n+1}$ and $p^{n+1/2}$ by solving a time-centered system of equations:

[TABLE]

where $\mathbf{f}^{*}$ is calculated from $\mathbf{X}^{*}$. Finally, we update the positions in a time-centered manner to obtain $\mathbf{X}^{n+1}$:

[TABLE]

where $\mathbf{U}^{*}$ is computed in terms of $\mathbf{X}^{*}$. Although this scheme involves the first-order accurate preliminary guess $\mathbf{X}^{*}$, the approach is similar to the broader class of Runge-Kutta schemes and it can be shown to be second-order accurate (see for instance the related scheme in [10, 11]).

To actually solve for the velocity and pressure in (45), we take the discrete divergence of both sides of the equation and use (46) to eliminate the discretely divergence-free velocity field (we assume the initial velocity satisfies $\mathbf{D}\cdot\mathbf{u}^{0}=0$). This results in a Poisson problem for the pressure $p^{n+1/2}$, and once this problem has been solved (45) may be rewritten as an explicit equation for $\mathbf{u}^{n+1}$. We solve this linear system using geometric multigrid [29, 30].

3.3 Computing the height function

Given the sensitivity of the underlying equations to the height, e.g. the dependence in (33) of $V_{h}-V_{0}$ on $h^{3}$, the accurate computation of the height function is essential to this method. This additional computational cost compared to the standard immersed boundary method, in which computing the height is not necessary, is the price paid for resolving velocity gradients across thin fluid layers.

To compute the height $h_{i}$ above a point on the lower surface with index $i$ and position $\mathbf{X}_{0,i}$ (Figure 2), we approximate the upper surface by a piecewise cubic curve. On an arbitrary section of the upper surface with endpoints $\mathbf{X}_{h,j}$ and $\mathbf{X}_{h,j+1}$, the curve is given by $\mathbf{X}_{h,j}(s)=\mathbf{a}_{0}+\mathbf{a}_{1}s+\mathbf{a}_{2}s^{2}+\mathbf{a}_{3}s^{3}$ for $s\in[0,1]$. The four unknown constants \left\{\mathbf{a}_{k}\big{|}k=0,\dots,3\right\} are determined by the conditions

[TABLE]

where $\mathbf{T}_{h,j}$, the tangent vector to the upper surface at point $j$, is defined as the weighted average of the two edge vectors connected to vertex $\mathbf{X}_{h,j}$:

[TABLE]

Solving the system (48)–(51) results in the coefficients

[TABLE]

Defining components of the tangent vector by $\mathbf{T}=(T^{1},T^{2})$, the unit normal $\mathbf{N}$ is given by

[TABLE]

To calculate the height above $\mathbf{X}_{0,i}$, we find the pair $(s,h)$ and the index $j$ such that

[TABLE]

On a section $j$ of the piecewise cubic upper curve, we find $s$ by taking the dot product of both sides with $\mathbf{T}_{0,i}$ to obtain

[TABLE]

and the resulting nonlinear equation for $s$ is solved using Newton’s method. Having found $s$, the height $h$ is given by

[TABLE]

To take advantage of the quadratic convergence of Newton’s method, it is crucial to have a good initial guess. We use a piecewise linear (rather than piecewise cubic) approximation to the upper surface to obtain the initial guess $s_{0}$. This has the virtue of being easy to solve: finding $s_{0}$ such that $\mathbf{X}_{0,i}+h_{0}\mathbf{N}_{0,i}=\mathbf{X}_{h,j}+s_{0}(\mathbf{X}_{h,j+1}-\mathbf{X}_{h,j})$ yields

[TABLE]

The above procedure is done for each $j$ corresponding to a section of the upper surface. For some applications, there is a unique section that yields $s\in[0,1]$ and $h>0$. However, in some cases there are multiple sections with positive heights above $\mathbf{X}_{0,i}$ and $s$ in the appropriate range, such as when the upper surface is a closed curve. In that case, the height is defined to be the minimum of the eligible values. We note that a naïve implementation of this algorithm to compute heights would involve a nested loop over $i$ and $j$ and result in an $\mathcal{O}(N^{2})$ operation. However, looping over all $j$ can be avoided since the lubricated method is only applicable when two immersed boundaries are separated by a distance comparable to the grid spacing. Only nearby points need to be checked, and this allows for fast methods based on binning points into containers of size on the order of $\Delta x$ by $\Delta y$.

We apply the lubrication equations (5)–(7) at a point $\mathbf{X}_{0,i}$ together with its corresponding point $\mathbf{X}_{h,j}(s)$ on the upper surface. Since in general $\mathbf{X}_{h,j}(s)$ does not coincide with a discrete node, this involves defining a velocity $\mathbf{U}^{\text{IB}}_{h,j}(s)$ and Lagrangian force $\mathbf{F}_{h,j}(s)$ at the position $\mathbf{X}_{h,j}(s)$. We interpolate these quantities from the discrete points on the upper surface by the same cubic scheme used above to compute $\mathbf{X}_{h,j}(s)$. Note that although we use the lubrication corrections to solve for both $\mathbf{U}_{0,i}$ and $\mathbf{U}_{h,j}(s)$, only $\mathbf{U}_{0,i}$ is actually used. The velocities on the upper surface are calculated in an analogous but separate step.

We also need to define heights $h_{j}$ below each point on the upper surface in order to compute velocities there. The procedure used to do so is similar to the one just described for the lower surface, although since the height is defined relative to the lower surface normal, the equation that needs to be solved is slightly more complicated. We describe this procedure in detail in Appendix A.

3.4 Outer derivatives

To compute outer derivatives of velocity, we use the global velocity field $\mathbf{u}^{\text{IB}}(\mathbf{X},t)$ and the fact that the discrete delta function $\delta_{h}(\mathbf{x})$ is continuously differentiable. The global velocity gradient is defined by

[TABLE]

where $\nabla\delta_{h}(\mathbf{x})=\left(\frac{1}{\Delta x}\phi^{\prime}(x/\Delta x)\phi(y/\Delta y),\frac{1}{\Delta y}\phi(x/\Delta x)\phi^{\prime}(y/\Delta y)\right)$ and $\phi^{\prime}(r)$ is calculated by analytically differentiating the expressions in (43) to obtain

[TABLE]

As shown in the numerical experiments of Section 4.1, the interpolated gradient gives a good approximation to the outer velocity derivative when it is evaluated at least 2 gridpoints from the boundary. We therefore calculate the outer derivative at a distance of 2 – 4 gridpoints from the immersed boundaries. The optimal choice depends on the application; for instance, a smaller distance is sometimes needed to avoid evaluating $\nabla\mathbf{u}^{\text{IB}}$ outside the domain.

4 Convergence Results

4.1 Simple shear

As a first demonstration of the lubricated immersed boundary method, we simulate the situation described in Section 2.2 of two parallel lines separated by a distance $h$ and pulled in opposite directions.

Figure 3: Test problem of two parallel lines pulled in opposite directions by tether points.

See Figure 4.1 for an illustration. We implement this scenario by using tether points, i.e. we prescribe the motion of virtual lines $\mathbf{Z}_{0}(t)$ and $\mathbf{Z}_{h}(t)$ and connect immersed boundaries to these lines by stiff springs with spring constant $k_{\text{tether}}$. Tether points are widely used in immersed boundary simulations to prescribe either rigid boundaries, as in Wiens et al. [31], or moving boundaries with prescribed kinematics, as in the simulations of jellyfish in Hamlet et al. [32]. The Lagrangian forces are given by

[TABLE]

In our simulations we use a periodic domain of size $L_{x}=L_{y}=2$, with a prescribed shear rate of $\dot{\gamma}=0.15$ and a height $h=1/24$ between the two lines. The fluid has a density of $\rho=1$ and a dynamic viscosity of $\mu=0.02$, which together with the characteristic velocity of $0.004$ gives a Reynolds number of $\text{Re}\approx 0.2$. The units are taken to be dimensionless. The timestep is refined according to $\Delta t=0.01\Delta x$. Each line is parameterized in the Lagrangian domain $[0,2\pi]$ and discretized using $N_{r}$ points so that $\Delta q=2\pi/N_{r}$. We choose $N_{r}=2N_{x}$ to satisfy the rule of thumb that, to prevent leakage through boundaries, the spacing between points on immersed boundaries should be approximately half the grid spacing [12]. Further, we increase the tethering constant together with the grid resolution according to $k_{\text{tether}}=12.5\cdot(N_{r}/16)^{2}$, so that the lines become rigid in the limit $\Delta x\to 0$. Simulations are run up to a total time of $T=10$.

Figure 4: Results from applying the standard immersed boundary method to the problem of parallel lines pulled in opposite directions. The asterisks in LABEL:sub@figtest1_ib:vel show the velocity values used to advect the lines. Note in LABEL:sub@figtest1_ib:l1err that all errors in $v$ and $p$ are $\lesssim 10^{-12}$, so that rounding errors dominate.

Given the tether point kinematics prescribed in our case, as $k_{\text{tether}}\to\infty$ we expect to recover $\mathbf{X}_{0}=\mathbf{Z}_{0}$ and $\mathbf{X}_{h}=\mathbf{Z}_{h}$ with a linear shear flow throughout the domain. The results obtained using the standard immersed boundary method are illustrated in Figure 4.1, with the mean velocities of the immersed boundaries denoted by asterisks. The standard method yields accurate results provided that $h\gg\Delta y$. However, as would be reasonably expected it becomes less accurate when $h\lesssim\Delta y$. Because of the underresolved relative velocity between the two lines, extreme spring forces are generated that lead to large overshoots in the fluid velocities near boundaries. See Table 4.1 for the errors obtained on different grid sizes, including the corresponding relative error in the mean shear rate $\dot{\gamma}=(\overline{U}_{h}-\overline{U}_{0})/\overline{h}$, where $\overline{U}_{h}$ and $\overline{U}_{0}$ are the average velocities on each surface and $\overline{h}$ is the average height.

Table 1: Convergence results for the standard immersed boundary method applied to the problem of parallel lines pulled in opposite directions. We give absolute $\textrm{L}^{1}$ and $\textrm{L}^{\infty}$ errors in the velocity components $(u,v)$, pressure $p$, and immersed boundary position $\mathbf{X}$ and velocity $\mathbf{U}$. We also report the relative error in the mean shear rate $\dot{\gamma}$.

The lubricated immersed boundary method, on the other hand, resolves the fluid flow through the thin gap exactly in the special case of linear shear and therefore converges rapidly to the correct velocities. As shown in Figure 4.1, the resulting fluid velocities do not exhibit the large overshoots that appear in the standard immersed boundary method. The errors computed by the lubricated immersed boundary method are significantly smaller. For instance, on a coarse $16\times 16$ grid for which $\Delta y=3h$, the errors made by the lubricated method are over an order of magnitude smaller (see Table 4.1 for $\textrm{L}^{1}$ and $\textrm{L}^{\infty}$ errors on different grid sizes and Movies flat_ib64.gif and flat_lube64.gif for a visual comparison of the two methods). Because lubrication corrections are applied in this case, the asterisks denoting the mean velocities of the immersed boundaries do not coincide with the interpolated velocities from the fluid grid.

Figure 5: Results from applying the lubricated immersed boundary method to the problem of parallel lines pulled in opposite directions. The asterisks in LABEL:sub@figtest1_lube:vel show the velocity values used to advect the lines. Because the lubricated immersed boundary method is used, the asterisks need not coincide with the interpolated velocities. Note in LABEL:sub@figtest1_lube:l1err that all errors in $v$ and $p$ are $\lesssim 10^{-13}$, so that rounding errors dominate.

On coarse grids that underresolve the fluid, using the lubricated immersed boundary method results in accuracies that would require a $2-4$ times finer grid if the standard immersed boundary method were used. Since the number of degrees of freedom scales with the inverse grid spacing squared, this implies that the lubricated method can obtain the same accuracy approximately $8-64$ times faster, taking into account both the number of degrees of freedom and the CFL stability restriction $\Delta t\lesssim\Delta x$.

Note that on a $128\times 128$ grid for which $\Delta y=3h/8$, the standard immersed boundary method gives nearly the same error as the lubricated method. This implies that when the gap contains more than 2 gridpoints, the linear shear flow is sufficiently resolved and the lubricated method is no longer necessary.

Table 2: Convergence results for the lubricated immersed boundary method applied to the problem of parallel lines pulled in opposite directions.

4.2 Eccentric rotating cylinders

The test problem of two parallel lines is a relatively simple case, of course. The assumptions of lubrication theory are actually valid for linear shear flow independent of the gap size, and the absence of curvature hides approximations made by the lubricated immersed boundary method. Therefore, we next consider the more demanding problem of eccentric rotating cylinders.

This problem is appealing because it is simple to describe and implement but gives rise to non-trivial fluid dynamics. Unlike the case of linear shear flow, the lubrication approximation is not exact in this case. Moreover, interesting phenomena emerge such as lift on the inner cylinder and counter-rotation of the fluid on the wide side of the gap, as we describe next. In discussing the setup of this well-known problem we closely follow the exposition of Acheson [24].

On a domain of size $L_{x}=L_{y}=2$, we use inner and outer cylinders of radii $r_{1}=3/4$ and $r_{2}=3/4(1+1/24)$ respectively, so that the nondimensional thickness $\epsilon:=(r_{2}-r_{1})/r_{1}=1/24$. The outer cylinder is centered at the origin, whereas the inner cylinder is centered at $(x_{0},0)$ with $x_{0}=3/128$ so that the cylinders are not concentric. With this choice of parameters, the dimensionless eccentricity $\lambda:=x_{0}/(r_{1}\epsilon)$ satisfies $\lambda=3/4$. This value is chosen partly because whenever $\lambda>(\sqrt{13}-3)/2\approx 0.30$, the exact solution exhibits a counter-rotating vortex on the wide side of the gap [24]. Moreover, there is a lift $F$ on the cylinder equal to $F=12\pi\mu U\lambda/(\epsilon^{2}\sqrt{1-\lambda}(2+\lambda^{2}))$.

Next, we investigate which of these features can be recovered on coarse grids using the lubricated immersed boundary method. To simulate this problem, we again use tether points with specified kinematics. The outer cylinder is taken to be stationary, while the inner cylinder is tethered to points rotating counterclockwise with prescribed angular velocity $U=8.33\cdot 10^{-4}$. The simulation is run up to a total time of $T=100$. We use a fourfold stronger spring constant of $k_{\text{tether}}=50\cdot(N_{r}/16)^{2}$ compared to the linear shear simulations of Section 4.1. All other parameters are unchanged.

Figure 6: Standard immersed boundary method applied to the problem of eccentric rotating cylinders. Note the absence of reverse flow in LABEL:sub@figtest2_ib:ugap, i.e. $u_{\theta}\geq 0$ everywhere. The gap in LABEL:sub@figtest2_ib:ugap has been blown up for purposes of visualization by scaling all distances from the inner cylinder by a factor of $1/\epsilon=24$.

On the coarsest $16\times 16$ grid used, the two cylinders come to within a minimal distance of $\Delta x/16$, and attain a maximal distance of about $\Delta x/2$. Given these subgrid separations, it is no great surprise based on the results of the previous section that using the standard immersed boundary method leads to large overshoots near the boundaries that are significantly improved by the lubricated method (Figures 4.2 and 4.2). Further, by computing the tangential velocity through the gap using (34), we are able on a coarse, $16\times 16$ grid to detect reverse flow indicating the existence of a counter-rotating vortex, whereas using the standard method we cannot observe this even using a significantly finer, $128\times 128$ grid. See Figures 4.2LABEL:sub@figtest2_ib:ugap and 4.2LABEL:sub@figtest2_lube:ugap for a comparison of the two methods on a $32\times 32$ grid.

The lubricated method accurately resolves the pressure through the gap (Figure 4.2LABEL:sub@figtest2_lube:press), allowing for the accurate calculation of lift. The lift is calculated by summing the Lagrangian force density over the inner immersed boundary and multiplying by $\Delta q$ to obtain a force. The lift values obtained using the lubricated immersed boundary method on $16^{2}$, $32^{2}$, $64^{2}$, and $128^{2}$ grids have corresponding relative errors 67%, 27%, 2.1%, and 0.58%, consistent with 2nd-order accurate convergence. The lifts calculated using the standard immersed boundary method, on the other hand, actually have the wrong sign for grids smaller than $128^{2}$, and converge to the correct value so slowly that there is still a $14\%$ error on a $512^{2}$ grid. These significant errors are not entirely surprising given that the standard method does not accurately capture discontinuities at the fluid-structure interface [33].

Interestingly, the advantages of the lubricated method are not made readily apparent by the errors in certain other quantities calculated at several grid sizes (Tables 4.2 and 4.2). The standard immersed boundary method is able to separate the cylinders nearly as well as the lubricated method and does not seem to exhibit a significantly underresolved relative velocity between the cylinders. Indeed, the mean shear rate is the only quantity displayed in Tables 4.2 and 4.2 that results in an error at least an order of magnitude smaller for the lubricated method. However, as in the linear shear test problem, there are large overshoots in the fluid velocities (see Movies ecc_ib64.gif and ecc_lube64.gif) and other quantities such as the gap pressure are not computed accurately. This shows the importance for assessing accuracy of investigating auxiliary quantities, such as the lift and flow reversal, that are not controlled directly through the tether points.

Table 3: Convergence results for the standard immersed boundary method applied to the problem of eccentric rotating cylinders. We give absolute $\textrm{L}^{1}$ and $\textrm{L}^{\infty}$ errors in the velocity components $(u,v)$, pressure gradient $\mathrm{d}p/\mathrm{d}\theta$, and immersed boundary position $\mathbf{X}$ and velocity $\mathbf{U}$, along with the relative error in the mean shear rate $\dot{\gamma}$ computed in terms of the angular velocities. The pressure gradient is computed by applying the jump condition (22) to obtain $p$ and then taking a centered finite difference along the boundary.

Figure 7: Lubricated immersed boundary method applied to the problem of eccentric rotating cylinders. Note the presence of reverse flow in LABEL:sub@figtest2_lube:ugap, i.e. $u_{\theta}<0$ on the far side of the gap. The gap in LABEL:sub@figtest2_lube:ugap has been blown up for purposes of visualization by scaling all distances from the inner cylinder by a factor of $1/\epsilon=24$.

Table 4: Convergence results for the lubricated immersed boundary method applied to the problem of eccentric rotating cylinders.

5 Applications

5.1 Wall-induced migration

Many fluid-structure interaction problems in biology, including blood flow, intracellular trafficking, and leaf aerodynamics involve the near-contact of multiple objects. Blood is a non-Newtonian fluid that exhibits several counter-intuitive behaviors, and one interesting such phenomenon is the emergence of a cell-free layer at the vessel walls in which no red blood cells are present. The cell-free layer has its origins in wall-induced migration, which causes deformable objects such as liquid droplets and red blood cells to acquire lift and move away from walls (see Figure 5.1LABEL:sub@fig_onewall:exp and [34, 35, 36, 37]). Because blood in the microcirculation is at low Reynolds numbers ($\text{Re}\approx 10^{-4}-10^{-2}$) so that the Stokes flow approximation is valid, directed motion away from the wall may seem to be at odds with the well-known scallop theorem that asserts no net displacement may be achieved by a time-symmetric motion in Stokes flow [38]. The paradox is resolved in this case since symmetry is broken as the object deforms with the flow [35, 36].

Figure 8: LABEL:sub@fig_onewall:exp Experimental demonstration of wall-induced migration of an elastic vesicle, reprinted from [37] with the permission of AIP Publishing, LABEL:sub@fig_onewall:sim_lube Wall-induced migration simulated by the lubricated immersed boundary method on a $64\times 64$ grid. The position of the elastic ring is shown at several times as it moves through the channel and away from the wall. LABEL:sub@fig_onewall:sim_ib Analogous simulation by the standard immersed boundary method; the ring fails to migrate from the wall because of the underresolved relative velocity between the two boundaries. Axes are given in the dimensionless units $\hat{x}=x/R_{0}$, $\hat{y}=y/R_{0}$.

Simulating wall-induced migration by the standard immersed boundary method is challenging because of the issues discussed thusfar in resolving the lubricating layer that separates the deformable object from the wall. However, by applying the lubricated immersed boundary method, motion from the wall can be seen even on coarse grids. We have studied this phenomenon by simulating an elastic ring moving through a channel. The elastic ring is made up of immersed boundary points $\mathbf{X}_{i}$ for $i=0,\,1,\,\dots,\,N_{q}-1$ arranged in a loop and connected by springs having total energy

[TABLE]

and by torsional springs with total bending energy

[TABLE]

where the indices are computed modulo $N_{q}$. We simulate the motion away from the wall using the standard and lubricated immersed boundary methods and compare the qualitative results (see Figure 5.1). The elastic ring has radius $R_{0}=\pi/4$, $k_{\text{spring}}=1$, and $k_{\text{bend}}=0.004$, and it is immersed in a periodic domain of unit cell size $L_{x}=4\pi$ by $L_{y}=2\pi$ that is filled with fluid having density $\rho=1$ and viscosity $\mu=0.2$. The viscosity is taken to be spatially uniform. This is the appropriate assumption for systems such as lipid-bilayer endosomes and red cell ghosts in which the internal and external fluids share the same material properties.

The channel walls are built by putting down a single line of tether points at $y=0$. One wall of tether points is sufficient to define the channel because the periodic boundary conditions in $y$ imply that a particle will reach the wall regardless of whether it moves up or down.

Flow is established in the channel by applying a body force of the form $(f_{x}=\mu\sin(y),f_{y}=0)$. In the absence of an elastic ring, this establishes a unidirectional flow $(u=\sin(y),v=0)$ in the domain. We use a $64\times 64$ grid (in which case $\Delta x\neq\Delta y$) and timestep $\Delta t=0.005\Delta x$, and 32 discrete points around the circumference of the ring. The elastic ring is initialized nearby the wall with its center of mass at $(x,y)=(0,R_{0}+\pi/192)$. This implies that at the beginning of our simulation the ring is at a distance of $\Delta y/6$ to the wall. The simulations are run up to a time of $T=28$.

In these simulations of elastic vesicles, we use a modified fluid solver that has been used previously to simulate inkjet printers [39] and is based on the finite element projection of [28]. A marker-and-cell (MAC) fluid velocity is constructed at each time step in which the horizontal and vertical velocities are defined at cell edges and are not co-located. The advantage of this fluid solver is that the MAC velocity is useful as an interpolating field since advected elastic structures tend to suffer less leakage and conserve volume better than those advected on collocated grids [40].

Figure 5.1 shows that, at a qualitative level, the lubricated immersed boundary method is able to capture the sequence of deformations during wall-induced migration. The elastic ring simulated by the lubricated immersed boundary method deforms and moves away from the wall, whereas the ring simulated by the standard immersed boundary method becomes stuck and cannot travel with the flow (see Movies onewall_ib64.gif and onewall_lube64.gif).

5.2 Channel flow

Figure 9: Comparison between simulations of channel flow using the standard (left column) and lubricated (right column) immersed boundary methods, with the color map representing the pressure. An elastic ring is initialized at a distance of $\Delta y/6$ from both walls and a uniform force is applied to the ring. In the standard immersed boundary method simulation, the ring sticks to the wall near its starting position, whereas in the lubricated immersed boundary method simulation the ring separates from the channel walls and moves in the direction of the applied force.

Figure 10: Comparison between simulations of channel flow using the standard and lubricated immersed boundary methods, with the color map representing horizontal velocities at time $t=40$. There is an order of magnitude difference in the horizontal velocities obtained; in the standard immersed boundary method simulation, the channel is effectively clogged.

Finally, we perform a simulation of an elastic ring flowing through a narrow channel. This example is inspired by the fluid-structure interaction that takes place in the microcirculation when red blood cells squeeze through capillaries as small as a few µm in diameter.

The size of the elastic ring in the following simulation is chosen so that it fits just within a channel defined by two walls of tether points. We use an elastic ring of radius $R_{0}=9\pi^{2}/40$ that is initialized at a distance of $\pi/192$ from each wall. As in the previous section, we use a $64\times 64$ grid on a domain of size $L_{x}=4\pi$ by $L_{y}=2\pi$ so that, in terms of the grid spacing, the ring is initialized at a minimal distance of $\Delta y/6$ from each wall. The ring is forced in the $+x$-direction by applying a uniform force to each immersed boundary point. Here we discretize the ring with $64$ immersed boundary points around its perimeter and use the timestep $\Delta t=0.01\Delta x$. The simulation is run up to $T=40$. Note that for this problem, we must compute two height functions, one for the distance of the vesicle to each wall.

When the lubricated method is used, the ring flows along in the channel, adopting bullet-like shapes reminiscent of red blood cells in capillaries. This is in contrast to results obtained by the standard immersed boundary method, in which the ring becomes stuck at its starting point (Figure 5.2). The horizontal velocity field computed by the standard method displays grid artifacts in the velocities in the narrow gap separating the elastic ring from the walls, and small velocities throughout the rest of the channel, indicating a clog. The horizontal velocities computed by the lubricated immersed boundary method, on the other hand, do not display these grid artifacts in the lubrication layers, and are over an order of magnitude greater in the rest of the channel (Figure 5.2). See Movies channel_ib64_p.gif and channel_lube64_p.gif, which show the evolution of the pressure field, and Movies channel_ib64_u.gif and channel_lube64_u.gif, which show the evolution of the horizontal velocity, for a comparison of the results obtained by the standard and lubricated methods.

6 Conclusions

We have proposed an immersed boundary method that uses aspects of lubrication theory to resolve thin fluid layers between elastic structures. The two key ingredients of the method are (i) the observation that across a lubrication layer the normal derivatives of velocity at the boundaries determine the relative velocity, and (ii) the well-established jump conditions that relate these normal derivative of velocity to the Lagrangian force on the boundary. We have applied the lubricated immersed boundary method to problems of increasing complexity and have found that the method outperforms the standard immersed boundary method when the lubrication layer is small compared to the fluid grid spacing. The lubricated method makes it possible on coarse fluid grids to capture interesting phenomena, such as the lift and reverse flow in the problem of eccentric rotating cylinders, and wall-induced migration in the problem of an elastic vesicle near a wall in shear flow.

Note that, although in the problems considered here the immersed boundaries were initialized to be close to touching, immersed boundaries that are initially well-separated can also come within subgrid distances over the course of a simulation. This effect is particularly apparent in underresolved 3D computations with significant differences in velocities from gridpoint to gridpoint, as occurs for instance in simulations of red cells moving through capillaries [10] and flagellar bundling in E. Coli [18]. As noted in the Introduction, in this last case initially separated immersed boundaries have even been observed to cross, making it a good candidate to which to apply the lubricated method.

Although we have demonstrated the effectiveness of the present formulation in several settings, much work remains to be done to make the lubricated immersed boundary method readily applicable to real-world problems. The results presented here are restricted to 2D, and although much of the mathematical formulation extends to 3D, it is not clear to what extent the specific numerical approach (e.g. the definition of the height function in terms of piecewise cubic curves) will extend to the challenges of 3D problems. Further investigation may indicate that alternative formulations are preferable to the one presented here; for example, the current formulation is not symmetric with respect to upper and lower surfaces, a property that would be desirable for some applications. (However, for those applications in which one surface tends to be significantly flatter than the other, an asymmetric formulation may be preferable.) Moreover, in many applications it may be more convenient to impose physical boundary conditions within the fluid solver rather than using a tethered surface. In future work, we plan to find appropriate lubrication corrections in these types of situations as well.

7 Acknowledgments

We thank Charles Peskin for useful conversations and the anonymous reviewers for their helpful feedback. TGF acknowledges funding through National Science Foundation grant DMS-1502851. CHR was supported by the Applied Mathematics Program of the U.S. Department of Energy (DOE) Office of Advanced Scientific Computing Research under contract DE-AC02-05CH11231.

Appendix A Computing the height on the upper surface

Consider a point $\mathbf{X}_{h,j}$ on the upper surface (Figure 11). We approximate the lower surface by a piecewise cubic curve as before. That is, given a section of the lower surface with endpoints $\mathbf{X}_{0,i}$ and $\mathbf{X}_{0,i+1}$, the curve between these endpoints is $\mathbf{X}_{0,i}(s)=\mathbf{a}_{0}+\mathbf{a}_{1}s+\mathbf{a}_{2}s^{2}+\mathbf{a}_{3}s^{3}$ for $s\in[0,1]$, where the coefficients are defined in an analogous manner to (48)–(51).

We wish to find the pair $(s,h)$ and the index $i$ such that

[TABLE]

where as before $\mathbf{N}_{0,i}(s)$ is defined by $\mathbf{N}=(-T^{2},T^{1})/\|\mathbf{T}\|$ through the tangent vector, which in this case is given by differentiating the cubic surface, i.e.

[TABLE]

Since the normal vector itself depends on the unknown $s$, the formulation (108) is subtly different from the analogous formulation (58) of the height above a point on the lower surface. To solve for $s$, we take the dot product of (108) with $\mathbf{T}_{0,i}(s)$ to obtain the 5th order polynomial equation

[TABLE]

and once more use Newton’s method to compute the solution. Given $s$, it is straightforward to solve for the height via $h_{j}=\left(\mathbf{X}_{h,j}-\mathbf{X}_{0,i}(s)\right)\cdot\mathbf{N}_{0,i}(s)$. As before, it is important to have a good initial guess to ensure the efficient convergence of Newton’s method. To obtain this initial guess, we approximate the lower surface by a piecewise linear function so that $\mathbf{X}_{0,i}(s)=\mathbf{X}_{0,i}+s\left(\mathbf{X}_{0,i+1}-\mathbf{X}_{0,i}\right)$, but we use the continuously varying normal vector

[TABLE]

The continuously-varying normal is required to ensure there is a region of finite width around the lower surface above which heights can be computed. Now, we wish to find $s_{0}$ such that

[TABLE]

Taking dot products with $\mathbf{T}_{0,i}(s_{0})$ yields the quadratic equation

[TABLE]

where

[TABLE]

Of the two roots $(-b\pm\sqrt{b^{2}-4ac})/(2a)$, the negative root is the relevant one since it remains finite in the limit of zero curvature in which $a\to 0$. We have found that using the quadratic formula to actually compute the roots can lead to a loss of precision; instead we use the numerically robust method described in [41].