Background on immersed boundary method

Algorithm for the Immersed Boundary Method: In an immersed boundary simulation, the state of the system at time $t_{n}=n\Delta t$ is given by the fluid velocity and pressure at that time; the locations of the immersed boundary points, $\mbox{$\bf X$}_{i,p}(t)$ on each immersed boundary entity i; and the configuration of elastic links which join these points. To update the system from time t_n to time t_n+1, the following basic steps are carried out:

1.

Calculate an elastic force $\mbox{$\bf F$}_{i,p}$ at each immersed boundary point $\mbox{$\bf X$}_{i,p}$.

2.

Transmit the immersed boundary forces $\mbox{$\bf F$}_{i,p}$ to the fluid, by spreading each $\mbox{$\bf F$}_{i,p}$ from the corresponding immersed boundary point $\mbox{$\bf X$}_{i,p}$ to the uniform grid used to solve the Navier Stokes equations. The cumulative contributions of all the immersed boundary forces $\mbox{$\bf F$}_{i,p}$ gives the grid force density $\mbox{$\bf f$}_{jl}$ at each grid point (x_j,y_l):

 

$$\mbox{$\bf f$}_{jl} ~=~ \sum_{i,p} ~ \mbox{$\bf F$}_{i,p} \u... ...mbox{$\bf x$}_{jl} ~-~ \mbox{$\bf X$}_{i,p})}_{{\rm Spreading}}$$ (1)

3.

Solve discrete versions of the Navier-Stokes equations using $\mbox{$\bf f$}_{jl}$as the driving force

to obtain the grid velocity field $\mbox{$\bf u$}_{jl}$.

4.

Move each immersed boundary point at a velocity obtained by interpolating the grid velocity to the immersed boundary point:

 

$$\frac{\mbox{$\bf X$}_{i,p}^{n+1} - \mbox{$\bf X$}_{i,p}^{n}}... ...{$\bf x$}_{jl} ~-~ \mbox{$\bf X$}_{i,p})}_{{\rm Interpolation}}$$ (2)

5.

If appropriate, create new elastic links and/or destroy some existing elastic links.

6.

If appropriate, modify the properties of some of the elastic links.

Each of these steps is described in some detail in the following sections. To pick appropriate parameters for a simulation, it is necessary that the next section is clearly understood. In the above, many of the occurrences of $\mbox{$\bf X$}_{i,p}$ lack an indication of time level. For the current version of IBIS, the force calculations and movement of immersed boundary points are done using explicit time discretization so each unlabeled $\mbox{$\bf X$}_{i,p}$ should be interpreted as $\mbox{$\bf X$}_{i,p}^{n}$, the values at the beginning ot the timestep. In the future, we will provide an option for doing these calculations in an implicit manner, in which case some or all of the unlabeled $\mbox{$\bf X$}_{i,p}$might be interpreted as $\mbox{$\bf X$}_{i,p}^{n+1}$, the unknown values at the end of the timestep.

Immersed Boundary Elastic Forces: To specify the elastic properties of the immersed boundaries, an energy function of the immersed boundary point coordinates is created for each entity (or set of connected immersed boundary points), and the elastic forces that act on the immersed boundary points are calculated from the derivatives of the energy function. For each immersed boundary entity, the energy function has three contributions; energy from stretching springs, energy from bending the entity, and energy from moving the entity away from a specified location. One or more of these contributions may be zero for a given entity at a given time.

We assume that the elastic links which generate the stretching energy act as Hooke's Law Springs. This means that if elastic link l connects immersed boundary points with indices p1(l) and p2(l) and with coordinates, $\mbox{\boldmath$X$}_{p1(l)}$ and $\mbox{\boldmath$X$}_{p2(l)}$, respectively, then the appropriate energy function for this link is:

 

$$E_S(\mbox{\boldmath$X$}_{p1(l)},\mbox{\boldmath$X$}_{p2(l)})... ..._{p1(l)} - \mbox{\boldmath$X$}_{p2(l)}\vert\vert - R_{l} ] ^2.$$ (3)

where R_l is the link's resting length, S_l is its stiffness coefficient and $\Vert \mbox{$\bf X$}\Vert$ denotes the length of the vector $\mbox{$\bf X$}$. This function is zero when $\Vert \mbox{\boldmath$X$}_{p1(l)} - \mbox{\boldmath$X$}_{p2(l)} \Vert = R_{l}$ and attempting to minimize this energy drives the distance between the immersed boundary points towards the resting distance.

It is important that the points of each immersed boundary entity be spaced closely enough so that no fluid leaks across the immersed boundary. For this it suffices that consecutive points on an immersed boundary be separated by a distance of no more than h/2, where h is the grid spacing. So, for links between consecutive points of the same entity, it is important that R_l < h/2. For elastic links between points of distinct entities, this constraint on R_l need not hold.

The bending energy involves a pair of links emanating from the same immersed boundary point, and is designed to penalize deviations in the angle between these links from a prescribed angle. Suppose that immersed boundary points $\mbox{\boldmath$X$}_{p}$ and $\mbox{\boldmath$X$}_{q}$ are joined by an elastic link, as are points $\mbox{\boldmath$X$}_{q}$ and $\mbox{\boldmath$X$}_{r}$.We define the bending energy of this triple of points to be:

 

$$E_B(\mbox{\boldmath$X$}_{p},\mbox{\boldmath$X$}_{q},\mbox{\b... ...imes (\mbox{\boldmath$X$}_{q}-\mbox{\boldmath$X$}_{p}) - C ]^2.$$ (4)

Here S_B is the bending stiffness, C is a specified `curvature', and $\hat{z}=(0,0,1)$. This function achieves its minimum value when the exterior angle $\theta$ between the links satisfies

 

$$\theta = \sin^{-1} \frac{C}{\Vert \mbox{\boldmath$X$}_{r} -\... ...Vert \Vert\mbox{\boldmath$X$}_{q}-\mbox{\boldmath$X$}_{p}\Vert}$$ (5)

If the angle between the pair of links differs from $\theta$,the force associated with the energy drives the angle toward $\theta$.If the resting lengths for the links in the triple are R_pq and R_qr, then in order to set the desired angle to $\theta_{\rm{desired}}$, C should be set to:

$$C = R_{pq} R_{qr} \sin(\theta_{\rm{desired}}).$$

Typically, the three points in the bending energy formula are consecutive points of a given immersed boundary entity, and one has a contribution of this kind for each triple of consecutive points along the boundary. The bending stiffness S_b and desired curvature C can vary from triple to triple.

The energies E_S and E_B are used to define the elastic properties of the immersed boundary entity. Terms like E_S help determine the size of the entity and prevent the immersed boundary points from separating far enough to allow flow to cross the immersed boundary. Terms like E_B help determine the shape of the immersed boundary; in particular, they are useful to prevent the boundary from becoming too irregularly shaped near any point.

The displacement or `tether' energy is designed to penalize deviations from a prescribed location. If immersed boundary point $\mbox{\boldmath$X$}_{p}$ is connected by a spring to tether point $\mbox{\boldmath$X$}_{p}^{T}$,the corresponding tether energy is:

 

$$E_T(\mbox{\boldmath$X$}_{p}) = \frac {1}{2} S_{T,p} \Vert \mbox{\boldmath$X$}_{p}- \mbox{\boldmath$X$}_{p}^{(T)} \Vert^2.$$ (6)

Here, S_T,p is the tether spring stiffness. The tether point $\mbox{\boldmath$X$}_{p}^{(T)}$ is not an immersed boundary point, and its motion can be prescribed. Often it is a fixed point in space, in which case the force derived from the tether energy drives the immersed boundary point towards this fixed point. Used in this way, the tether points and tether energy can help fix the absolute location of the immersed boundary points in the face of an externally-driven flow. In some applications, tether points can be used to cause an immersed boundary to move in an approximately prescribed manner. For example, if immersed boundary point $\mbox{\boldmath$X$}_{p}$is connected to a tether point $\mbox{\boldmath$X$}_{p}^{T}$ by a link with large stiffness coefficient S_T,p, and the tether point is moved, say, to the right at speed s, then the immersed boundary point will tend to move right at a speed of approximately s.

The elastic energy of entity i is the sum of the stretching, bending and tether energies for each immersed boundary point in that entity. Suppose, for example, that entity i consists of a string of N immersed boundary points arranged in the order $\mbox{\boldmath$X$}_{i,1},\mbox{\boldmath$X$}_{i,2},\ldots,\mbox{\boldmath$X$}_{i,N}$; that each pair of consectutive points is joined by an elastic link; that each consecutive triple has an associated bending energy; and that each point $\mbox{\boldmath$X$}_{i,p}$ is linked to a tether point $\mbox{\boldmath$X$}^{T}_{i,p}$. Then the elastic energy of this entity is:

 

$$E_i(\mbox{\boldmath$X$}_{i,1},\mbox{\boldmath$X$}_{i,2},\ldo... ...$}_{i,q+1}) \,+\, \sum_{q=1}^{N} E_T(\mbox{\boldmath$X$}_{i,q})$$ (7)

Here, $\delta = 0$ if there is a spring connecting $\mbox{\boldmath$X$}_{1,i}$with $\mbox{\boldmath$X$}_{N,i}$, and $\delta = 1$ otherwise. In the former case, which corresponds to a closed string, we treat the indices of the immersed boundary points as periodic with period N, and so identify points with indices and N, 1 and N+1, and so on. Other examples might also include stretching energy from elastic links which connect immersed boundary points in entity i to immersed boundary points in a different entity, as well as bending energy from triples of linked points on two entities.

The total elastic energy for the system at time t is the sum of the elastic energies of the individual entities of the system:

 

$$E(\mbox{\boldmath$X$},t) ~=~ \sum_{i=1}^{M} E_i(\mbox{\boldm... ...\mbox{\boldmath$X$}_{i,2},\ldots,\mbox{\boldmath$X$}_{i,N_{i}})$$ (8)

where M is the number of entities in the system, N_i is the number of immersed boundary points in entity i, and $\mbox{\boldmath$X$}$ is short-hand for the set of coordinates of all of the immersed boundary points. The elastic force at point p of entity i, $\mbox{\boldmath$F$}_{i,p}$ is then calculated from the derivatives of the elastic energy $E(\mbox{\boldmath$X$},t)$with respect to the coordinates in $\mbox{\boldmath$X$}_{i,p}$:

 

$$\mbox{\boldmath$F$}_{i,p}(\mbox{\boldmath$X$},t) ~=~ - {\par... ...E(\mbox{\boldmath$X$},t)}/{\partial \mbox{\boldmath$X$}_{i,p}}.$$ (9)

The minus sign is chosen to drive the system towards a minimum of the energy.

Notice that to specify the energy and the forces associated with an entity, the constants S_l, R_l, S_B, C, and S_T,p must be specified. This can be accomplished in several ways, as is discussed below. We recall also that these quantities can vary from entity to entity and from point to point, or link to link, within an entity. How the user specifies this is also discussed below.

Spreading the forces to the finite difference grid: We now suppose that all of the immersed boundary forces $\mbox{\boldmath$F$}_{i,p}$have been computed for this timestep, and we discuss how these are transmitted to the grid used for the solution of the Navier-Stokes equations. Since the immersed boundary points need not coincide with the grid points, we require a means of communicating information between these two sets of points. To this end, we define the function $D_{jl}(\mbox{$\bf X$})$ where j and l are indices of a point of the fluid grid, and $\mbox{$\bf X$}=(X,Y)$ is an arbitrary point in space:

 

$$D_{jl} (X,Y) ~=~ \delta_h (X - x_j) ~ \delta_h (Y - y_l)$$ (10)

where

$$\delta_h (x) = \left\{ \begin{array} {ll} \frac{1}{4h} ~(1 ... ...{ } & \mbox{ } \ 0 & \mbox{, otherwise} \end{array} \right.$$

The function $D_{jl}(\mbox{$\bf X$})$ has the property that

 

$$\sum_{j,l} D_{jl}(\mbox{$\bf X$}) h^{2} ~=~ 1$$ (11)

for any point $\mbox{$\bf X$}$. Note also that $D_{jl}(\mbox{$\bf X$})$ is nonzero only in a $4h \times 4h$ portion of the grid surrounding $\mbox{$\bf X$}$.The fluid force density is defined from the immersed boundary forces using this function $D_{jl}(\mbox{$\bf X$})$ as follows:

 

$$\mbox{$\bf f$}^{n}_{jl} ~=~ \sum_{i,p} \mbox{$\bf F$}^{n}_{i,p} D_{jl}(\mbox{$\bf X$}^{n}_{i,p})$$ (12)

In this sum, the function $D_{jl}(\mbox{$\bf X$})$ is evaluated at each of the immersed boundary point locations $\mbox{$\bf X$}^{n}_{i,p}$ in turn. Because of the property (12) of $D_{jl}(\mbox{\boldmath$X$})$, the sum of $\mbox{$\bf f$}^{n}_{jl}$ over all points of the grid equals the sum of $\mbox{\boldmath$F$}^{n}_{i,p}$ over all immersed boundary points, all of the immersed boundary force affects the fluid motion. Also, because $D_{jl}(\mbox{$\bf X$}_{i,p})$ is zero except for a small region centered on $\mbox{$\bf X$}_{i,p}$, the immersed boundary forces are spread by (13) to a thin layer of the domain along each immersed boundary.

Solution of the Navier-Stokes Equations: The computational solution of the Navier-Stokes equations is currently done using a finite-difference method introduced by Chorin. Chorin's scheme is an example of what is known as a `projection' method. In general, projection methods take into account fluid advection, viscous forces, and outside forces (here the forces included in $\mbox{$\bf f$}$) to define a provisional velocity field which does not necessarily satisfy the incompressibility condition $\nabla \cdot \mbox{$\bf u$}= 0$. Then this provisional velocity is projected on the set of incompressible velocity fields, giving both the final velocity at the end of the timestep and the pressure. Chorin's method was one of the first projection methods, and has many nice features. It also has some limitations, one of which effectively restricts our calculations to relatively small Reynolds numbers (e.g., less than 50 or so) at the present time. We are working to implement a more recent projection method without this limitation. When this is ready, we should be able to handle Reynolds numbers of up to about 1000.

Moving the immersed boundary points: To move the immersed boundary points according to (3) we have to interpolate the fluid velocities from the fluid grid to the location of each immersed boundary point. To do this we again make use of the function $D_{jl}(\mbox{\boldmath$X$})$ defined above. Let $\mbox{$\bf U$}^{n}_{i,p}$ be defined by

$$\mbox{$\bf U$}^{n}_{i,p} ~=~ \sum_{jl} \mbox{$\bf u$}^{n}_{jl} D_{jl}(\mbox{\boldmath$X$}^{n}_{i,p}) h^{2}.$$

Thus, Uⁿ_i,p is an average of the velocity values at fluid grid points in a $4h \times 4h$ region around the location $\mbox{\boldmath$X$}^{n}_{i,p}$. The location of immersed boundary point is then updated using the formula:

 

$$\mbox{\boldmath$X$}^{n+1}_{i,p} ~=~ \mbox{\boldmath$X$}^{n}_{i,p} ~+~ \Delta t U^{n}_{i,p}.$$ (13)

This is done for each immersed boundary point.

Dynamic Link Creation and Removal: In some situations, such as modeling cell-cell cohesion or adhesion of a cell to a substrate one would like to allow new elastic links between entities to develop during the course of a simulation, and to allow these same links to break under certain conditions. Two such situations which have been studied are platelet aggregation during blood clotting, and biofilm deposition. What one needs in these situations is a set of rules for deciding when links should be formed and when they should be broken, and code that facilitates tracking these dynamic links.

To make our description concrete, we will consider the case of Fogelson's platelet simulations. In these simulations, immersed boundary strings are used to model the two walls of a blood vessel and other nearly circular immersed boundary entities are used to model the platelets suspended in the blood. The walls and platelets can interact and links between them can be formed or broken.

For example, whenever two `activated' platelets are sufficiently close together, an elastic link may be formed to connect a point on one platelet ring to a point on the other platelet ring. This simulates cohesion. At most one such link may emanate from any particular point on a platelet, and the maximum number of links allowed between a given pair of platelets is specified. Whether, in a given timestep, a new link is formed between two platelets depends on the distance between them and the current number of links connecting them. So, tests are done between pairs of (nearby) eligible platelet points to see if they are close enough to be linked. Once formed, a cohesive link whose index is l makes a contribution of the form

 

$$E_{{\rm coh}(l)} ~=~ \frac{1}{2} \, S_C \left[ \Vert \mbox{\... ...),p_1(l)} - \mbox{\boldmath$X$}_{i_2(l),p_2(l)} \Vert \right]^2$$ (14)

to the energy function. Here, we assume that this link connects point p₁(l) on platelet i₁(l) to point p₂(l) on platelet i₂(l). The spatial coordinates of these points are $\mbox{\boldmath$X$}_{i1_l,p1_l}$ and $\mbox{\boldmath$X$}_{i_2(l),p_2(l)}$, respectively. The resting length of each cohesion link is assumed to be . Similar links are formed when a platelet ring adheres to the injured part of the vessel wall. In the case of adhesive links, the energy contribution is

 

$$E_{{\rm adh}(l)} ~=~ \frac{1}{2} \, S_A \left[ \Vert \mbox{\... ...,p_1(l)} - \mbox{\boldmath$X$}_{i_2(l),p_2(l)} \Vert \right]^2.$$ (15)

Here i₁(l) is the index of the platelet and p₁(l) is the point on that platelet from which the adhesive link emanates, and i₂(l) and p₂(l) are the indices of the wall and the point on the wall, respectively, to which the platelet adheres.

Cohesive and adhesive links are also permitted to break. Testing for this involves looping over all cohesive and adhesive links and determining whether the separation between the points joined by the a given link exceeds a prescribed yield-distance. In this case, the link is removed from the calculation. The sets of cohesive links L_c(t) and adhesive links L_a(t) are dynamic. They grow if new links are formed, and shrink if links are broken, and these are tracked using linked-list-based data structures in the platelet aggregation calculations. Recall that the we assume that the fluid motion is spatially periodic with period the size of the computational domain. By default the immersed boundaries are also treated in a periodic fashion, so that if its motion carries an immersed boundary entity out of, say, the right end of the domian, it reappears at the corresponding location at the left end of the domain. For some situations this is not desirable. For example, in the platelet aggregation simulations, the properties of a platelet may change due to its contact with the wall or exposure to fluid-borne chemicals. When this happens the platelet is said to be `activated'. Normally this happens because of some local injury to a blood vessel and so platelets approaching the injury from upstream are not activated, and some of those which move downstream of the injury have become activated. It would be very unrealistic in this context, to allow the activated platelets that move downstream to reenter the domain at the upstream end. So for these calculations, we use a routine to remove platelets that reach the downstream end of the computational domain. We also use a routine to randomly introduce new platelets near the upstream end of the domain at a prescribed rate. A consequence of this is that the number of immersed boundary entities in the platelet calculations may change over the course of the simulation.

We have not yet generalized sufficiently our procedures for dynamic link formation and breaking, and for dynamic immersed boundary creation and removal to include them in the current version of IBIS, but we plan to do so in the near future.

Actively changing the elastic link properties: In some applications, it is important that the elastic properties of the immersed boundary change with time. For example, in modeling muscle contraction using a set of elastic links, the resting lengths of the links would shorten at a prescribed speed during the active contraction of the muscle. In modeling swimming by aquatic flagellates, propagation of a wave of bending along a flagellum can be accomplished by appropriately varying the curvatures C in (5) at different immersed boundary points along the flagellum. Many other such applications can be imagined. IBIS provides a routine for extracting immersed boundary point and spring data and passing it to a user-supplied subroutine udefine in which the user implements rules for dynamically changing the spring properties and/or tether point locations.