Immersed Interface Methods for Stokes Flow with Elastic Boundaries or Surface Tension.
A second-order accurate interface tracking method for the solution of incompressible Stokes flow problems with moving interfaces on a uniform Cartesian grid is presented. The interface may consist of an elastic boundary immersed in the fluid or an interface between two different fluids. The interface is represented by a cubic spline along which the singularly supported elastic or surface tension force can be computed.   The Stokes equations are then discretized using the second-order accurate finite difference methods for elliptic equations with singular sources developed in our previous paper [SIAM J. Numer. Anal., 31(1994), pp. 1019--1044]. The resulting velocities are interpolated to the interface to determine the motion of the interface.  An implicit quasi-Newton method is developed that allows reasonable time steps to be used.
Introduction
. In this paper we develop an interface tracking method for the solution
of incompressible Stokes flow problems with moving interfaces on a uniform Cartesian
grid. The interface may consist of an elastic boundary immersed in the fluid, as in the
model problem of Tu and Peskin[49], or an interface between two different fluids, as in the
study of bubbles or free surfaces.
The method we use is based on the Immersed Interface Method (IIM) for elliptic problems
developed in our previous paper [22] and the second author's thesis[26]. This is a
second order accurate Cartesian grid method for solving elliptic equations whose solutions
are not smooth across some interface, due to discontinuous coefficients or singular source
terms in the equation. The main idea is to incorporate the known jumps in the solution or
its derivatives into the finite difference scheme, obtaining a modified scheme whose solution
is second order accurate at all points on the uniform grid even for quite arbitrary interfaces.
This approach has also been applied to parabolic equations[28], [30] and hyperbolic wave
equations with discontinuous coefficients[23], [24].
Similar ideas have been used in the context of domain embedding, where an irregular
region is embedded into a larger rectangular domain on which a Cartesian grid is used.
Methods of this type include the method of fictitious domains[2], [8], [16], [29], capacitance
matrix methods[9], [43], and Mayo's method[32], [33].
A variety of Cartesian grid methods have been proposed for fluid dynamics problems
with arbitrary boundaries and/or moving interfaces, e.g., [1], [4], [11], [12], [17], [21], [31],
[37], [44], [45], [51], [50], [52] and other references below. Such methods are becoming
increasingly popular for problems with very complex geometries or moving interfaces where
more standard "body-fitted" or unstructured grid approaches can run into difficulties.
This work was supported in part by NSF Grants DMS-8657319, DMS-9204329, DMS-9303404, DOE
Grant DE-FG06-93ER25181 and URI Grant N00014092-J-1890.
y Departments of Mathematics and Applied Mathematics, University of Washington, Seattle, WA 98195.
(rjl@amath.washington.edu).
z Department of Mathematics, UCLA, Los Angeles, CA (zhilin@math.ucla.edu).
Our motivation for the present work was the goal of ultimately developing a second
order accurate version of Peskin's Immersed Boundary Method, a very robust algorithm
for solving the full incompressible Navier-Stokes equations with moving boundaries. This
method was originally developed for studying blood flow in a beating heart[38], [39], [40],
but has also been used in a wide variety of other problems (e.g., [6], [13], [14], [51]) particularly
in biomechanics since these problems often involve complex geometries and are
difficult to model with more standard approaches. The physical domain is embedded in
a rectangular region and the Navier-Stokes equations are solved on a uniform Cartesian
grid. The boundary is modeled by a singular forcing term which, computationally, is then
approximated by a set of discrete delta functions which spread the force from the boundary
to the nearby Cartesian grid points. This is described further in the context of Stokes flow
in [49]. See [39] for complete details. The spreading of forces via discrete delta functions
appears to limit Peskin's method to first order accuracy. Our approach will, we hope, lead
to second order accurate methods for the Navier-Stokes equations.
In this paper we consider the simpler Stokes equations which are still of physical interest
in many applications. The Stokes equations model the creeping flow of a highly viscous fluid,
in the limit where the Reynolds number goes to zero and both the inertial acceleration and
convection terms are dropped from the Navier-Stokes equations. We will concentrate on
describing the algorithms in two space dimensions. In Section 11 we briefly discuss what is
involved in extending them to three space dimensions.
In two dimensions the Stokes equations take the form
(1.1a)
(1.1c)
Here ~u j (u; v) is the velocity vector, p is the pressure, - is the viscosity, and ~
is the external force. We also use the notation below. See, for example, [3], [10],
[19], [20], [42] for general discussions of Stokes flow.
The equations (1.1) can be solved as a coupled system (as is done in [49]) or alternatively
reduced to a sequence of three Poisson problems, one for each variable. Differentiating (1.1a)
with respect to x, (1.1b) with respect to y, and adding the equations together gives
F
where r 2 is the Laplacian. Since the right hand side is known, this is a Poisson problem
for the pressure. Once p is known, the equations (1.1a) and (1.1b) are independent Poisson
problems for u and v respectively.
Note that the time evolution of the flow is governed entirely by the time-dependence
of the forces ~
F . If ~
F is known at a given instant in time, then the system (1.1) is elliptic
and the solution (u; v; p) is determined independently of the past history of the flow. This
is a reflection of the fact that there is no inertia in the system, i.e., the convective and
time-derivative acceleration terms have been dropped from the momentum equations. This
allows us to use the techniques developed in [22] directly, once we have determined the
appropriate jump conditions on the solutions. We assume that the reader is familiar with
that paper.
The jumps in the solution result from the fact that the force ~
F is singular and is
supported only along the interface (resulting from the elastic force or surface tension). This
singular force in (1.1a) and (1.1b) leads to jumps in the first derivatives of u and v across
the interface. Since the Poisson problem (1.2) for p involves derivatives of ~
F , and hence a
dipole source, the pressure will be discontinuous, along with its derivatives.
The singular force can be written as
~
Z L 0~
X(s;
where ~
X(s; t) gives the location of the interface at time t, parameterized by s for
~
f(s; t) is the force strength at this point, and ffi is the two-dimensional delta function. This
singular force is best viewed as a distribution whose action on any smooth test function
OE(~x) is
Z L 0~
X(s; ds:
In practice we eliminate the singular source term from the right hand side of (1.2) and
instead solve
together with specified jump conditions across the interface, using the techniques of [22].
This is described in more detail in Section 4.
We first consider the model problem used by Tu and Peskin[49], an immersed elastic
band in a fluid with the same fluid properties on each side. In Section 10 we extend our
approach to study an interface between two different fluids, with surface tension providing
the force at the interface rather than an elastic band.
The elastic band problem is the two-dimensional analog of an elastic balloon in a highly
viscous fluid (with the same fluid inside and out). In equilibrium, an ideal balloon would
have a spherical shape, with zero velocity everywhere and uniform pressure both inside
and out, but with a jump in pressure across the elastic membrane. The magnitude of this
jump depends on how far the membrane is stretched from its resting configuration. In two
space dimensions the analog is an elastic band, or more physically the cross section of an
impermeable cylindrical elastic tube, which contains an incompressible fluid and is stretched
to a diameter greater than its resting diameter. The equilibrium configuration is a circle
with a jump in pressure that balances the elastic force exerted by the stretched membrane
in a manner made clear in Section 3.
Note that we make an important distinction between the equilibrium configuration and
the resting configuration. By the equilibrium configuration we always mean the steady
state solution with a given quantity of fluid inside the elastic band, which will typically be
stretched like an inflated balloon since the fluid cannot escape. By the resting configuration
we mean the shape that would ultimately be obtained if a small leak were introduced,
allowing the band to deflate until there is no force exerted by the elastic membrane and no
pressure difference across it. The length of the elastic band in this unstretched state will
be called the resting length, L 0 , and the radius of the corresponding circle will be denoted
by r inflated band, at equilibrium, will have some equilibrium length L e and
radius r depend solely on how much fluid is trapped inside.
If we perturb an inflated band from its equilibrium position, analogous to squeezing
an inflated balloon, it will return to an equilibrium circular shape. Our goal is to simulate
this motion. We specify an initial position of the boundary (some simple closed curve) and
simulate its motion to a circle. The area A enclosed by the band must be invariant with
time and the equilibrium radius must be r
A=-.
We will parametrize the location of the band at any time t by ~
where s is arclength along the unstretched band, measured from some arbitrary
but fixed origin. The material particle at (X(s; t); Y (s; t)) is the same particle that would
lie a distance s from the origin if the band were cut at the origin and allowed to relax. The
force exerted by the band at ~
X(s; t) is given by (1.3) with strength
~
@s
where T (s; t) is the tension at this point, given by
@ ~
@s
and ~- (s; t) is the tangent vector to the band at this point,
@s
@ ~
@s
Note that in the relaxed state j@ ~
since s is the arclength. The single scalar T 0
describes the elastic properties of the band, which are assumed to be uniform along its
length, though this is not necessary. The larger the stiffer the elastic band and the
greater the force induced by a stretching of the band. The tension given by (1.7) is a linear
Hooke's law relation, which could easily be replaced by a more general nonlinear relation
in our algorithm.
The force density can be computed directly from the location of the interface ~
X at time
t. We will see in Section 3 that it is possible to determine the jump conditions for u; v and
p directly from ~
f(s; t). This in turn allows us to apply the immersed interface method to
solve for u; v and p at all points on a uniform Cartesian grid at time t.
The interface at each time is represented by a cubic spline passing through a given set
of control points along the interface, as described in the next section, and hence the location
of the interface is completely determined by the location of the control points. Taking a
time step requires also moving the interface, which is accomplished by moving the control
points using the additional constraint that the interface must move with the fluid. The
velocity ~u will be continuous across the interface and we have the differential equation
@
@t
~
The velocities computed on the Cartesian grid are interpolated to the control points, which
are then moved with this velocity over a time step of length \Deltat. In practice an implicit
method must be used in order to take reasonable time steps, and a quasi-Newton method
has been developed to accomplish this in an efficient manner, as described in Section 7.
This completes an overview of our method. Each of these steps will be described in
greater detail in the remainder of the paper.
Our approach differs from that of Tu and Peskin in several ways. The fundamental
difference is our use of jump conditions to derive second order accurate finite difference
schemes, rather than using discrete delta functions (which appear to limit the accuracy to
first order). But there are other differences as well.
Using a relatively small number of control points to specify the interface rather than
the denser set typically used with the Immersed Boundary Method has great advantages
in solving the implicit system of equations. This approach could be carried over to the
Immersed Boundary Method too, as described in Section 8.
We also use a sequence of three Poisson problems as described above, whereas Tu and
Peskin discretize the Stokes equations (1.1) directly as a coupled system, which is necessary
in using the Immersed Boundary Method since the right hand side of (1.1) involves delta
function forces, whereas the right hand side of (1.2) has dipoles, which are even more difficult
to discretize accurately. Our approach, on the other hand, could also be applied to the
coupled system directly, incorporating the jump conditions for u; v, and p simultaneously.
This has also been implemented and details are presented in [26]. We have found, however,
that the approach presented here gives slightly better results. It also yields the pressure p
at all grid points, which may be of interest, whereas discretizing the coupled system and
solving by fast Fourier transforms, as done in [49], does not yield p due to zero divisors in the
transformed equations. Finally, the approach based on three Poisson problems allows the
direct use of software already developed for the general Poisson problem with jumps across
an interface, a problem that also arises in many other contexts. However, decoupling the
problem into three Poisson problems depends on the fact that periodic boundary conditions
are used. Some comments are made in Section 11 on handling other boundary conditions.
Another approach to solving Stokes flow problems is to solve a biharmonic equation
for the stream function instead of the three Poisson problems used here. Mayo[31] has
developed a method based on this approach and techniques from [32] that is similar in
spirit to our method.
Boundary integral methods are also very popular for Stokes flow, since this linear problem
can be reduced to an integral equation along the interface, reducing the dimensionality
of the problem. For a description of this approach and some references, see e.g., [42]. This
approach does not appear to extend to nonlinear problems such as the full Navier-Stokes
equations, however, whereas the Immersed Boundary Method does. We expect that our
method can also be extended, and work is now underway to do so.
2. Representation of the interface and forces. The location of the interface at
time t n is represented by a finite set of control points (X n
Since the boundary is always a simple closed curve in the model problem, we assume
). The kth control point gives an approximation to (X(s
. Based on these control points, we determine a continuous
curve (X n (s); Y n (s)) by computing some interpolants X n (s) through the points X n
and Y n (s) through the points Y n
k .
The computations presented here were all computed using periodic cubic splines for
(s). The force exerted by the elastic interface, given by (1.6), can also be
represented by cubic splines. From the splines (X n (s); Y n (s)) we can compute @ ~
X=@s and
hence the tension T (s; t n ). Multiplying by the tangent vector and differentiating again gives
an approximation to (1.6). We now evaluate this function at each of the control points to
obtain values ~
k and then interpolate these values by a new cubic spline ~
f n (s) to obtain
the representation of the force all along the interface. The reason for introducing a new
cubic spline at this point is that the jump conditions discussed in the next section depend
on further derivatives of the force along the interface.
Other representations of the interface are possible. We have also experimented with
using a Fourier series representation, which is quite convenient for the smooth closed curves
considered here since X and Y are both periodic in s. With this approach it is easy
to compute the necessary derivatives and also apply filtering to remove high frequency
oscillations of the interface. In some problems this has been found to improve stability
properties.
A level set representation (e.g., [36], [47]) could also be considered. In this approach
the interface is represented by the OE(x; contour of some smooth function OE that
advects with the fluid velocity. For the surface tension problems discussed in Section 10,
this could be used since the surface tension force depends only on the curvature of the
interface, which can be determined directly from OE. However, for the elastic band problem
the force depends on the manner in which the band is stretched and not just its location,
and hence it seems necessary to explicitly track control points on the interface.
3. Jump conditions across the interface. In order to use the Immersed Interface
Method, we need to know the jump conditions for each of the three Poisson problems (for
u, and v). We need to know both the jump in the function and the jump in its normal
derivative at each point along the interface. (See [22] for details on how the jump conditions
are used to derive second order accurate difference formulas.)
The velocity components u and v are both continuous across the interface. Certainly
the normal velocity must be continuous and a discontinuity in the tangential velocity is also
ruled out by the presence of viscosity and a no-slip boundary condition between the elastic
band and the fluid on each side. The normal derivative of all three variables will, however,
be discontinuous in general, as will the pressure itself.
The jump conditions are easiest to write in terms of normal and tangential components
of the force ~
f(s; t). We decompose the force as
~
where ' is the angle between the x-axis and the normal direction pointing outward from the
interface at ~
X(s;
f 2 are the normal and tangential force strengths, respectively,
defined as
The jump conditions are then given by
(3.12a)
@s
(3.12c)
(3.12d)
These are derived in the Appendix.
Note that in the equilibrium case we expect there to be no tangential component to
the force, i.e., -
since in equilibrium the tension T (s; t) will be constant in s and so
~
f(s; t) will be the derivative of a constant-length tangent vector, and hence will point in the
normal direction. The only nonzero jump will then be that given by (3.12a), as we expect
since is piecewise constant in this case.
The jump conditions can be determined at any point on the interface by applying (3.12)
to the cubic spline function ~
representing the force. Alternatively, we could calculate
the jump conditions at each control point and then interpolate these by cubic splines.
4. Solving the Poisson problem for p. Given the location of the interface at time
and the jump conditions for p across it from (3.12a) and (3.12b), we wish to solve the
Poisson problem (1.5) with these specified jump conditions. We solve a difference equation
of the formh 2
on a uniform Cartesian grid, where the right hand side C ij is determined using the techniques
of [22] based on the jump conditions. The value C ij will be nonzero only at "irregular" grid
points, those for which the 5-point stencil includes points from both sides of the interface.
These values can be thought of as giving a spreading of the dipole source in (1.2) from the
interface to the nearby grid points, in the spirit of Peskin's Immersed Boundary Method.
However, they are determined from the jump conditions rather than the dipole in such a
way that the computed values p ij are second order accurate at all points (provided only
that the interface is smooth).
For the model problem of Tu and Peskin, we have periodic boundary conditions, and
so we know that a solution exists only if the C ij satisfy
0:
Unfortunately, since the method of [22] is based on requiring only that the local truncation
error be O(h) at irregular points, the values of C ij will have O(h) errors which, when
summed over the O(1=h) irregular points, will give an O(1) result for
C. However,
if we perturb each C ij to
we obtain a solvable perturbed systemh 2
whose solution is the least squares solution to the unperturbed equation (see [48]). Notice
that -
C=N 2 is of order h 2 . This means the order of the local truncation errors, which are
away from the interface and O(h) near the interface, have not been changed. We
expect the solution to the perturbed equation (4.14) will approximate the true solution of
the Poisson equation to second order accuracy. We use the Fourier method described in
[48] to solve the perturbed equation (4.14) and set p to get a particular solution.
5. Solving for u and v. Next we need to solve the Poisson problems (1.1a) and (1.1b)
for u and v. These are essentially identical and we will explain the procedure for the solution
of (1.1a) to obtain u. The forcing term F 1 is singular along the interface, but is a delta
function singularity rather than a dipole, leading to a jump in the normal derivative of u
but not in u itself. Also note that, since p is discontinuous across the interface, p x will
contain a delta function singularity. The strength of the jump in the normal derivative u n
depends on both these effects, and is given by (3.12c). As in our solution of the Poisson
problem for p, we do not attempt to model these singular terms directly, but instead drop
them from the equation and solve the Poisson problem
away from the interface, where p x is smooth, coupled with the jump conditions
(3.12c) across the interface. Using the IIM, this gives the discrete systemh 2
and the correction terms D ij , resulting from the jump conditions,
are nonzero only at the irregular grid points. These are again derived following [22]. Note
that if - is discontinuous across the interface then the elliptic equations for u and v have
discontinuous coefficients and the techniques of [22] can be used to derive a 6-point stencil
that yields second order accuracy. This is used in Section 10.
At regular grid points we can approximate p x by the standard central difference,
At irregular grid points we normally use one-sided differences to obtain first order accurate
approximations. Alternatively one could use the known jump conditions for p and p x to
correct the centered formula, but this also gives only a first order accurate formula unless
higher order jump conditions are derived and used. Fortunately, first order local accuracy
at these points is sufficient to maintain second order accuracy globally[22].
There is one situation where a one-sided approximation cannot be used, if both the
points lie on the opposite side of the interface from (i; j). This could
happen, for example, at the top or bottom of a circular interface. In this case we interpolate
to first order by using the known jump conditions [p], [p x ],
and [p y ] from (3.12a) and the relations
We know the jump [p n ] already from (3.12b), and
@s
@s
@s
i.e., [p s ] is the derivative of jump of pressure p along the interface. So we find
@s
@s
sin ';
@s
sin
@s
cos ':
be the control point closest to
the one of x i and x i+1 which is closer to X k . Then we can use the following interpolation
(5.
where the sign in the expression depends on which side of the interface the point (i;
on, and [p], [p x ] and [p y are calculated at (X k In [26] it is shown that (5.19) gives a
first order accurate approximation to p x .
Solving the system (5.15) gives velocities u ij that are second order accurate at all grid
points. An analogous procedure on (1.1b) gives the y-velocities v ij .
6. Moving the interface - an explicit method. The interface should move at
the local fluid velocity, which can be accomplished by moving each control point (X n
at velocity (U n
determined by interpolating the velocity fields u ij and v ij to the kth
control point. This interpolation is complicated by the fact that u and v are not smooth
across the interface, and so the known jumps in their normal derivatives must be used in
the interpolation formulas. There are various ways in which this can be done. Here we give
one example based on using 3 nearby points to perform linear interpolation, modified to
incorporate the jump conditions at the interface.
Dropping superscripts and subscripts for simplicity, let (X; Y ) be an arbitrary point
on the interface and consider the problem of interpolating from the u ij to obtain the x-component
U of the velocity at (X; Y ).
First we choose the first three grid points closest
to (X; Y ). Then we form a linear combination of the grid values at these points plus a
correction term to approximate U
We use Taylor expansion about (X; Y ) to get the equations for the coefficients fl's so that
we have a second order approximation:
Here the a i are linear combinations of the fl's as in equation (3.26) of [22], obtained from
the Taylor series expansion of each We then desire
This gives a linear system of equations for the fl's which can be solved to yield
Once we get the coefficients are able to compute the correction term C
as
We can use the same coefficients fl i and a new correction term based on the jumps in v to
find the y-component V of the velocity at (X; Y ). These velocities (U; V ) can then be used
to move the interface at the control point (X; Y ).
Applying this procedure at each control point (X n
k ) gives the velocities (U n
The simplest explicit method is forward Euler, in which we move the interface by shifting
each control point according to
Y n+1
This completes the description of an explicit immersed interface method for the Stokes
equations. In the next time step the whole process is repeated. To review, the process
consists of:
1. Use the location of the interface, as determined by the control points, to determine
the forces and jump conditions.
2. Solve three Poisson problems, using these jump conditions, to determine u n
ij and
ij on the uniform grid.
3. Interpolate to determine U n
k at the control points.
4. Move the control points at these velocities for time \Deltat.
There are two obvious difficulties. One is that Euler's method is only first order accurate
in time. A more serious difficulty is that the system is very stiff (for reasonable values of T 0 ),
and very small time steps must be taken to maintain stability, time steps over which there
is barely any discernable motion of the interface. The problem is that small perturbations
in the smoothness of the interface can lead to large forces, resulting in large transient
velocities that amplify the perturbations catastrophically if used over too large a time step.
This difficulty is discussed in detail by Tu and Peskin[49]. (See also [18], [34].) In order
to take reasonable time steps, we must use an implicit method, as described in the next
section.
7. Moving the interface - an implicit method. Steps 1 through 3 of the procedure
described in the previous section can be used to define an operator U that maps a set of
control points ~
to the resulting velocities ~
at the control points. We write
~
Applying U to ~
computing forces and jumps along the interface, solving three
Poisson problems, and interpolating the resulting velocities back to the control points. The
Forward Euler method of the previous section can now be written succinctly as
~
Instead we wish to use an implicit method, such as the Trapezoidal method
~
\Deltat (U( ~
This means that the distance each control point moves should be determined by the average
of its velocity based on the old interface location and its velocity based on the new interface
location. This is second order accurate and also eliminates most stability problems, but of
course it is much more difficult to apply. At time t n , ~
X n is known and so ~
can be computed as before. But then ~
must be determined from the implicit system
\Deltat U( ~
~
\Deltat ~
Normally a Newton-like method must be used with the Trapezoidal method in order to
obtain convergence of the iteration with reasonable size time steps. Newton's method
requires the Jacobian matrix
Unfortunately the matrix U 0 ( ~
X) is impossible to calculate exactly, and even obtaining a
finite difference approximation would be prohibitively expensive. To see this, recall that ~
represents the positions of all the control points (a vector with 2N b entries for N b control
points), while U( ~
X) represents the 2N b -vector of resulting velocities. Hence U 0 ( ~
X) is a
b \Theta 2N b matrix and its approximation by finite differences would require evaluating U
times. Since each evaluation of U requires solving three Poisson problems (plus other
work), this is out of the question.
Instead we use a quasi-Newton method, in which we maintain an approximation to the
Jacobian matrix J that is modified in each iteration by a low-rank update. The update
is based on information obtained in each step about the directional derivative of U in the
direction moved in this iteration. A wide variety of quasi-Newton methods are available,
see e.g., [15]. In practice one wants to avoid factoring the new matrix J in each iteration,
and so either the LU factorization is updated directly or else the inverse matrix
updated. We have implemented the BFGS (Broyden-Fletcher-Goldfarb-Shanno) method
with updates to B and have found this to be quite effective.
The efficiency of this method is enhanced by the fact that we have a very similar system
to solve in each time step. Since the interface moves only a small amount from one time
step to the next, the approximate Jacobian from one time step is still a good approximation
in the next time step. In each time step we begin with the final approximate B from the
previous time step as our initial B. Of course we also have a good initial guess for the
solution ~
X n+1 to the system g( ~
We can use ~
our initial guess, or even an
extrapolated value such as 2 ~
We find that we need only 2 or 3 iterations of the
quasi-Newton method to converge to ~
The above comments are valid once the process is going. In the very first time step,
from ~
X 0 to ~
we do not yet have a previous approximation to the Jacobian. In the first
step we initialize B to the identity matrix, which is reasonable since from (7.22) we see that
O(\Deltat). (This seems to work in spite of the fact that, in the stiff case, we want to
take \Deltat for which -\Deltat may be large for some eigenvalures - of U 0 .) In the first few time
steps, before a good approximate Jacobian has been built up, more iterations are required
than indicated above (5 or 6, typically).
8. Comparison with the Immersed Boundary Method. The original Immersed
Boundary Method can also be applied to Stokes flow, as has been done by Tu and Peskin[49].
We have already explained the primary differences between the two algorithms and, in
particular, how we expect to achieve second order accuracy rather than first order accuracy
by using jump conditions rather than discrete delta functions. These expectations are
confirmed in the numerical results presented in the next section.
There is another fundamental difference in approach that is worth discussing further
since we believe it could also be used to good effect in connection with the standard Immersed
Boundary Method. In Peskin's approach the Lagrangian points defining the boundary
are also the points where the discrete delta functions are centered in the process of
transfering the singular sources from the interface to the neighboring grid points. Since
the diameter of the support of these delta functions is O(h), where h is the uniform grid
spacing, this requires a fairly dense set of points along the interface. Typically the distance
between marker points on the interface is roughly h. With an implicit method, this can
give a very large system of equations to solve for the new position of the interface in each
time step.
By contrast, we use a much smaller set of control points to mark the interface, and
interpolate by cubic splines to determine intermediate locations or forces. If the interface is
smooth it is possible for the distance between these control points to be quite large relative
to h with no effect on the accuracy. This leads to a much smaller system of equations to
solve and allows us to use a fully implicit method quite efficiently.
Presumably improvements could be made to a code based on Peskin's approach by using
a hybrid method, in which the boundary is represented by a smaller set of control points
for the purpose of identifying its location (and hence in the implicit system to be solved).
One could then use cubic splines to interpolate forces to a denser set of points along the
interface where the discrete delta functions are applied for the purpose of spreading forces
and interpolating velocities. In essence this is what we do, since we also need to know the
forces and jump conditions at a denser set of points in the process of calculating the C ij in
the systems (4.13) and (5.15).
There are other advantages to solving for a smaller set of control points besides the
obvious efficiency gained by having a smaller system of equations. Reducing the number of
degrees of freedom in the interface position also leads to a much better conditioned system
in many cases. Beyer[5],[6] used the Immersed Boundary Method to model the motion of
the basilar membrane in the inner ear and implemented an implicit method to deal with
the rigid walls of the cochlea. He found that the system of equations was extremely ill-conditioned
and had to resort to a singular value decomposition to obtain good results, at
considerable expense. We suspect that using a much smaller set of control points would
eliminate the bulk of this ill-conditioning and lead to a much more efficient method for
handling rigid walls. Similar issues are also discussed in [18].
9. Numerical results. In this section we present the results from some sample calculations
on the Stokes equations with immersed elastic boundaries. We compare with results
obtained using the Immersed Boundary Method where the equations (1.1) are discretized
directly using a discrete delta function for the singular forces. This is the method used by
Tu and Peskin[49] and we have attempted to implement the method exactly as described in
that paper. For brevity we refer to that method as the IBM and to our Immersed Interface
Method as the IIM.
The main interest in [49] was in studying the issue of explicit versus implicit methods
for moving the interface in the context of a relatively simple test problem. While their
emphasis was not on obtaining optimal accuracy for this particular problem, we feel that
the comparison is appropriate since we hope to extend our method to the full incompressible
Navier-Stokes equations. It is important to verify that at least on this special case we are
able to obtain better results than the existing Immersed Boundary Method.
Example 9.1. This example is taken from Tu and Peskin [49]. The initial interface (the
solid line in Fig. 1) is an ellipse with major and minor axes a
The unstretched interface (the dash-dot line in Fig. 1) is a circle with radius r
to the restoring force, the ellipse will converge to an equilibrium circle (the dashed line
in Fig. 1) with radius r
ab - 0:61237; this is larger than the unstretched interface
because of the incompressibility of the enclosed fluid. So the interface is still stretched at
the equilibrium state, and the non-zero boundary force is balanced by a nonzero jump in
the pressure (Fig. 2(b)).
We begin by computing the velocities and pressure at time based on the initial
elliptical interface, before the interface has moved at all. Figure 2 shows u and p over the
uniform grid. As expected, p is discontinuous across the interface while u is continuous but
not smooth. Figure 3 shows this more clearly, displaying cross sections of u and p along
the line In Figure 3(b) we see that the discontinuity in pressure is captured very
sharply by the immersed interface approach.
In

Figure

3(a) we have also plotted the cross section of the velocity u that was computed
using our implementation of the IBM on the same grid. This gives a similar result except
near the interface, where it is quite smeared out with the sharp peak in velocity being lost.
Since it is the velocity right at the interface that is used to move the interface, this can be
expected to have a substantial impact on the overall performance of the algorithm. (We
cannot directly compare the accuracy of the pressure from the two algorithms, since this is
not available from Tu and Peskin's method.)

Table

1 shows the results of a grid refinement study on the IIM where the values on
three different N \Theta N grids with are compared to a fine grid solution
with The errors in p; u, and v are measured in the max-norm over all N 2 grid
points and displayed along with the ratios of successive errors. Since we are comparing to
a computed solution on a grid that is not very much finer, we do not expect the standard
error ratios of 2 for a first order method or 4 for a second order method. In particular, when
going from expect a qth order accurate method to produce a ratio
rather than the ratio 2 q . For ratio is 3 while for is 5. Table 1 shows the
final ratio to be betwen 5 and 6 for all three variables, indicating second order accuracy.

Table

shows results for the IBM, for the velocity components only. Now the final
ratios are roughly 3, indicating the expected first order accuracy.
We now consider the error at later times, after the interface has moved. Comparing the
solution at all the uniform grid points is difficult, since the interface may lie to one side of a
given point in one calculation, but slightly to the other side in a different calculation. Instead
we will focus on the interface location as a measure of the accuracy, which is appropriate
since this is often what we are most interested in.
One simple and effective measure is to study the values r min and r max , defined as
These measure the smallest and greatest distance from the origin to the interface. Note that
since we expect the interface to converge to a circle centered at the origin (by symmetry
and the positioning of the original ellipse), we expect that r n
-0.4
Initial
Final
Resting
Fig. 1. The interface at different states: The initial interface is the solid ellipse with
The equilibrium position is the dashed circle with re =
ab \Delta. The resting circle, shown as a
dash-dot line, has radius
(a)
-0.50.5-0.06
-0.020.020.06
x
y
(b)
x
y
pressure
Fig. 2. (a) The x-component of velocity u in the Stokes flow at It is continuous but not smooth
across the interface. (b) The computed pressure distribution of the Stokes flow at The pressure is
discontinuous.

Figure

4(a) shows how r min and r max behave computationally over a short time scale.
The solid line is from our IIM with and the dotted line is using the
IBM on the same grid. To plotting resolution these are indistinguishable on this scale.

Figure

4(b) shows what happens over a longer time scale, blown up near the true equilibrium
position r e - 0:61237. With our method, a numerical equilibrium is reached that agrees
well with the true equilibrium, and this equilibrium is then maintained. At 2000 we
have r min - 0:61232 and r max - 0:61248 and this is maintained at later times. With the
Immersed Boundary Method, some leaking is apparent which causes the circle to shrink.
This continues over longer times as well. At the circle has shrunk to a radius
r - 0:57 (not shown) and it still appears to be shrinking linearly with time. This problem is
also mentioned by Tu and Peskin[49]. Peskin and Printz[41] have suggested a modification
to fix the leaking problem, but we have not tried to implement this modification.
We can compare the accuracy over shorter times by again comparing the computed
solutions to a fine grid solution. The results for r max are shown in Figure 5 at two different
x
The x-component of the velocity u(x,-0.4).
(b)
-0.4
-0.20.20.611.4
x
Pressure profile: p(x,0)
Fig. 3. (a) A slice of u, the x-component of velocity in the Stokes flow at
is continuous but ux should be discontinuous across the interface. Solid line: IIM results, dot-dashed line:
IBM results. (b) A slice of pressure at The points and solid line both show the computed
results with the IIM at the grid points. (Pressure is not available with the IBM.) Note that the large jump
in pressure across the interface is computed with no smearing.

Table
The errors in the computed p, u, and v at using the IIM method via three Poisson equations.
Second order accuracy can be observed.

Table
The errors in the computed u and v at using the IBM. First order accuracy can be observed.
times, At the smaller time we see a clear improvement with our method.
At the later time our method still gives smoothly decreasing errors with a good rate. The
Immersed Boundary Method gives more chaotic results which look better on coarser grids
but poorer on fine grids. (The results seen are quite sensitive to the particular time chosen.)
Above we have considered the errors in r max . We could also look at some norm of the
error along the entire interface, for example the 2-norm. Let us take N
b as the finest
grid. For the coarser grid with N \Theta N , we take N In
this way we are guaranteed that each control point
x (N)
on the
interface is also a control point for the finest grid N   \Theta N   and N
b . Then it is possible to
(b)
20000.6020.6060.610.6140.618Fig. 4. Plot of rmax (upper curve) and rmin (lower curve) showing the distance from the origin to the
interface as a function of time t on a 160 \Theta 160 grid with on the boundary. Solid line: IIM results.
Dotted line: IBM results. (a) To plotting accuracy the behavior is identical over the time period
during which convergence to a near-circle is apparent. (b) Over a longer time scale, it is seen that the IBM
suffers leaking and the circle shrinks. Here the vertical scale is also expanded, showing that the numerical
equilibrium is not exactly circular in either case.
(a)
-5
log10(h)
o: IBM
(b)
-5
log10(h)
log10(E)
Fig. 5. Grid refinement studies of the error in rmax on a log-log scale. The solid line and the stars (
IIM results. The dash-dot line and 'o': IBM results. (a) At
compute the error
x (N)
i\Lambdal
y (N)
i\Lambdal
Our test results show that the error defined above is indeed a monotone decreasing function
as N increases. In Fig. 6, we plot the global error at with the finest grid being
being the pairs of (40; 40), (50; 40), (60; 40), (70; 40), (80; 80),
Note that the number of control points N b does not
decrease smoothly with N . We do this so that N b always divides 320, allowing direct
comparison at control points with the fine grid solution, as required in (9.23).
Fig. 6 shows that our method converges with a smaller error and a faster convergence
rate than the IBM. (Again the slope is greater than we would expect if we had the exact
solution to compare with, rather than a fine grid solution.)
-5
log10(h)
o: IBM
Fig. 6. The error in the interface location at measured in the norm (9.23). Solid line and the
star ( results. Dash-dot line and the 'o's: IBM results.
Another interesting phenomena we can observe from Fig. 6 is that the number of control
points N b on the interface plays an important role in the IBM. When we refine the uniform
grid but keep the number of control points N b on the interface fixed, the errors obtained
with their approach gradually cease to decline even as we refine the grid because the error in
expressing the interface will dominate. Then a refinement of the interface grid will lead to
a relatively large fall in the error as we can see in Fig. 6. There is a sharp drop in the error
with the IBM when N b changes from 40 to 80 and from 80 to 160. As expected, we must
refine the grid for the domain and the interface simultaneously with this approach. In our
approach, as we mentioned in Section 2, we can take fewer control points on the interface
with little effect on the accuracy with our method, as we can see from Fig. 6, where we
have the same sequence of grids. The jumps in the error when N b is increased are much
less pronounced.
Example 9.2. We now change the example slightly and perturb a circular interface
that is initially the resting circle, with r . There is still a restoring force that should
bring the interface back to a circular shape, but this force now vanishes as the equilibrium
(= resting) state is reached. Asymptotically there is no force exerted by the interface and
no jump in pressure. With both methods, the interface will stop moving before an exact
circle is reached, once the error in the discretization dominates the force. The deviation
from circularity of the final numerical equilibrium gives an indication of the accuracy of
the method. In Fig. 7 we have again plotted r min and r max , the minimum and maximum
distance from the origin to the interface, for both methods on a sequence of grids. It is clear
that our method is more accurate and it has been confirmed that r
large t, whereas with the Immersed Boundary Method, r
Example 9.3. This example shows that we can handle more complicated regions. The
initial interface in polar coordinates is 8'. The unstretched interface is the
(a)
x 100.550.650.75
time
(b)
x 100.550.650.75
time
Fig. 7. Plots of rmax (upper curves) and rmin (lower curves) for the case when the equilibrium state is
the same as the resting circle, re . Solid line, dash-dot line and dotted line are the results obtained with
160 \Theta 160, 80 \Theta 80 and 40 \Theta 40 grids respectively. (a) IIM results. (b) IBM results.
circle with the radius r present results only for 160. The
problem is very stiff and we need to take a fairly small time step even with the implicit
method. We started with \Deltat = O(h 2 ), but could increase the time step at later times. A
comparison with the IBM reveals a similar behavior as in Example 9.1. We will not give
detailed numerical results but instead present only the location of the interface at several
times in Fig. 8.
Fig. 8. The interface at different times with a 160 \Theta 160 grid. The dotted circle is the unstretched
interface with results only.)
10. Multi-fluid flows and surface tension. The method developed in the previous
section can be easily adapted for the interface between two different fluids, with surface
tension providing the singular force rather than an elastic membrane.
The force strength ~
f(s; t) is now given by
~
~
where the constant fl depends on physical properties of the two fluids, and s is arclength
along the interface. The vector @ 2 ~
X=@s 2 is normal to the interface with magnitude equal
to the curvature.
The main new feature that we need to include is the effect of gravity, which is important
in most applications since the two fluids may have different densities. If gravity is directed
in the negative y direction, we need only modify equation (1.1b) to read
where g is the gravitational constant and ae is the density, which we assume has the constant
value ae 1 in one fluid and ae 2 in the other. The Poisson problem for p then becomes
Since ae is piecewise constant, the term gae y gives only an additional delta function source
along the interface, which contributes to the jump in the normal derivative p n across the
interface. In fact this is the only contribution to this jump, since the force (10.24) is normal
to the interface and so -
(3.12). So we have
where the sin ' term arises from the fact that the delta function source is directed vertically,
and hence at angle ' to the interface. The jump conditions for p, u n , and v n are still given
by (3.12) with -
Note that the velocity is now continuously differentiable
across the interface, simplifying the procedure for interpolation to the interface that was
presented in Section 6.
Even if ae the addition of gravity will induce a hydrostatic pressure gradient
that is linear in y. This means that periodic boundary conditions are no longer reasonable.
However, if we are computing on the
and we set
Z Z
\Omega
ae(x; y) dx dy;
then we can write p as
where ~
p is the deviation from the linear profile obtained from the average density ae 0 . If the
boundaries are well away from the interface, then we expect ~
p to be roughly constant along
the entire boundary of @
\Omega\Gamma so the periodic boundary conditions are physically reasonable.
In terms of the pressure deviation ~
equation (10.25) becomes
~
Example 10.1. We consider a rising bubble of fluid with density ae inside the
bubble and density ae We take the same viscosity in both fluids, -
though this is not necessary (see Example 10.2 below). Figure 9 shows experiments with
different values of the surface tension coefficient and 0. In each case the
bubble was initialized to an elliptical shape,
When fl is large, the surface tension is sufficiently strong to bring the bubble back to a
nearly circular shape even as it rises. For smaller values of fl, the bubble is distorted. For
sufficiently small values, the bubble would eventually split into pieces. Our current code
cannot handle this change in topology. (See, e.g., [47] for such a calculation.)
This behavior agrees qualitatively with the known behavior of axisymmetric three-dimensional
bubbles, the case most frequently treated in the literature. (See for example
[3, 10, 46].) Our computations, however, are purely two-dimensional, which corresponds
to the cross-section of a cylindrical bubble in three dimensions. The additional curvature
effects in the more realistic axisymmetric case are known to have an effect on the shape of
bubbles (e.g., [35]) and so a direct comparison is not possible.
We also note that if we start with a circular bubble, rather than an ellipse, the bubble
remains circular (to reasonable accuracy) for all values of fl. This also agrees with expected
behavior[10].
Example 10.2. The next example shows that discontinuities in viscosity - can also
be incorporated. In this case the Poisson problems to be solved for the velocities u and v
have discontinuous coefficients and are solved using the techniques of [22]. The velocities are
continuous but now have jumps in their normal derivatives arising from the jump conditions

Figure

shows the results with inside, and for comparison also
the computation with As expected, the bubble moves more slowly in the
more viscous fluid. In this example we have also changed the boundary condition relative
to the previous example, and impose at the top and bottom, due to limitations in
our current code for the discontinuous coefficient Poisson problem. This explains why the
bubble with rises more slowly than the corresponding bubble in Figure 9.
These examples demonstrate that our technique is capable of dealing with discontinuities
in density and viscosity as well as singular forces at the interface. The ability to handle
discontinuous coefficients as well as singular forces is another advantage of our approach
over the standard Immersed Boundary Method, with which the discontinuities are typically
smeared over some region of width O(h), e.g., [47], [51].
11. Extensions. Although we have presented our method in the context of two-dimensional
flows, the ideas extend quite simply to three space dimensions. The Immersed
Interface Poisson solver for three-dimensional problems is described by Li[25],[26] and can
be used directly. The main difficulty is in representing the interface, which is now a surface,
by a finite set of control points. Various techniques are currently under consideration and
we hope to develop a three-dimensional code in the future.
The special case of axisymmetric flow can be solved as a two-dimensional problem, if
the Laplacian operator is suitably modified. In (r; y) coordinates, the Laplacian is r
r @ r +@ yy and this must be used in each of the three Poisson problems solved with our
approach. The jump conditions must also be modified to incorporate the three-dimensional
curvature effects.
Fig. 9. Bubble computations with different values of the surface tension parameter fl. In these computations
inside the bubble and everywhere. The computations are on a 160 \Theta 160
grid with points on the boundary. The computation with no surface tension
before when the bubble splits into pieces.
As noted in the introduction, decoupling the problem into three Poisson problems
depends on the fact that periodic boundary conditions are used. With other boundary
conditions, e.g., no-slip conditions at solid walls, we would know u and v but not p along
the boundary, and hence would not have a boundary condition for the Poisson problem
for pressure. If arbitrary boundary conditions are specified for p and the three Poisson
problems solved as described above, then the divergence will satisfy Laplace's equation
but this guarantees that u x +v y j 0 only if we also impose the boundary condition u x +v
along the boundary. This is the correct additional boundary condition which must be used
to uniquely define the pressure. This creates a troublesome coupling between the Poisson
problems.
One possibility is to introduce the values of p along the boundary as another set of
unknowns, with a discrete form of the divergence boundary condition as another set of
Fig. 10. Bubble computations with discontinuous viscosity. On the left
On the right, everywhere. The surface tension parameter cases.
equations to be solved via the quasi-Newton procedure along with the new locations of
the control points. This would give a much larger nonlinear system. Another possibility
is to revert to solving the coupled system of equations (1.1) directly as done by Tu and
Peskin. Results reported in [26] indicate that other aspects of the immersed interface
method presented here can be carried over directly to solving the coupled system.
In this paper we have considered only Stokes flow in the zero Reynolds number limit
of creeping flow, where the inertial time derivatives of the velocity can be neglected. We
are currently working on a version that includes these inertial terms in the momentum
equations (1.1a) and (1.1b). The techniques presented here can be used to compute the
time derivative of the velocities at both time t n and t n+1 , and then a Crank-Nicolson time
differencing can be used to achieve a second order accurate method in time. Again the
quasi-Newton method can be used to solve for the new interface location as part of this
updating process. The primary new difficulty arises from the fact that the time derivative
will not be smooth in time at any Cartesian grid point crossed by the interface during the
time step. This can be accounted for by including additional terms in the Crank-Nicolson
algorithm to correct for these jumps. Such techniques have been successfully used for one-dimensional
problems in [7], [27] and should carry over to multi-dimensional problems. This
work will be reported on elsewhere.
12.

Appendix

- Derivation of jump conditions. Here we present a brief derivation
of the jump conditions (3.12). We let ffi (~x) be the two-dimensional delta function as
defined in (1.3) and (1.4), and let OE(x; y) be an arbitrary twice continuously differentiable
test function defined on an appropriate
region\Omega\Gamma For a vector function ~
using Green's integral theorem, we
know that
Z Z
\Omega
Z
@\Omega
G \Delta ~n) OE ds \Gamma
Z Z
\Omega
~
We can use this to generalize the one dimensional result
Z
Fig. 11. Diagram used for the derivation of jump relations
to two dimensions, where r \Delta ~
F now behaves like the derivative of a delta function in the
direction normal to the interface. Multiplying by OE and integrating gives
Z Z
\Omega
Z Z
\Omega
aeZ
r \Delta ~
ds
oe
Z
@x
@y
ds:
Referring to Fig. 11, we take a belt
domain\Omega ffl which encloses the interface \Gamma.
ffl be the outer and inner boundary
of\Omega ffl respectively. Here ffl is some measure of the
distance from \Gamma to \Gamma
From the Poisson equation (1.2) for pressure we have
Z Z
Z Z
F dx dy
Z
@x
@y
ds:
Referring to Fig. 11, let us handle the left hand side of (12.27) first by using Green's theorem
repeatedly
Z Z
Z
ds
Z
ds \Gamma
Z Z
Z
ds \Gamma
Z
ds \Gamma
Z
ds
Z
ds
Z Z
where the superscripts indicate the values taken from outside and inside of the
interface \Gamma respectively. Notice that OE is twice continuously differentiable and p is bounded
and discontinuous only along the interface. So as ffl approaches zero, we have
Z Z
Z Z
Z
ds \Gamma
Z
ds as ffl ! 0:
To deal with the right hand side of (12.27), we express @OE=@x and @OE=@y in terms of normal
and tangential derivatives along the interface
@x
@y
sin ';
@x
sin
@y
cos ':
After solving the linear equations above for @OE=@x and @OE=@y, we get
@x
@y
Using these two equations in (12.27) and collecting terms, we have
Z
@x
@y
ds
Z
@n
@s
ds
Z
ds:
Integrating by parts in the second term and noting that \Gamma is closed, we may rewrite the
equation above as
Z
@x
@y
ds
Z
@s
OE
ds:
Comparing this with (12.28) and using the fact that OE is arbitrary, we must have
@s
To get the jump for u n , we multiply by OE(x; y) on both sides of (1.1a) and integrate
Z Z
Z Z
Z
ds:
The first term of the left hand side of (12.33) is
Z Z
Z
ds \Gamma
Z
ds \Gamma
Z Z
\Gamma!
Z
The second term of the left hand side of (12.33) is
Z Z
Z Z
p#
dx dy
Z
OE
ds
Z
OE
ds
Z Z
@x
@y
p#
dx dy
\Gamma!
Z
is arbitrary, from (3.11) and (12.35), we must have
Similarly, for v we can get
This completes the derivation of (3.12).

Acknowledgments

. We have benefited from conversations with many people during
the course of this work. In particular it is a pleasure to acknowledge the encouragement
and advice of Charlie Peskin. Much of this work was done while the second author was
a graduate student in Applied Mathematics at the University of Washington. This work
was completed while the first author was visiting the Scientific Computing Division of the
National Center for Atmospheric Research, supported by the National Science Foundation.



--R

A Cartesian grid projection method for the incompressible Euler equations in complex geometries.
Methods of fictitious domains for a second order elliptic equation with natural boundary conditions
An Introduction to Fluid Dynamics
A rotated difference scheme for Cartesian grids in complex geometries.
A computational model of the cochlea using the immersed boundary method
A computational model of the cochlea using the immersed boundary method
Analysis of a one-dimensional model for the immersed boundary method
On finite element domain imbedding methods
The direct solution of the discrete Poisson equation on irregular grids

Euler calculations for multieliment airfoils using Cartesian grids.
An adaptively-refined Cartesian mesh solver for the Euler equations
Interaction of oscillating filaments - a computational study
A mathematical model and numerical method for studying platelet adhesion and aggregation during blood clotting
New York
A fictitious domain method for Dirichlet problem and appli- cations
A fictitious domain method for external incompressible viscous flow modeled by Navier-Stokes equations
Removing the stiffness from interfacial flows with surface tension.
The Mathematical Theory of Viscous Incompressible Flow

High resolution finite Volume
The immersed interface method for elliptic equations with discontinuous coefficients and singular sources
Immersed interface methods for wave equations with discontinuous coef- ficients
Finite difference methods for wave equations with discontinuous coeffi- cients
A note on immersed interface methods for three dimensional elliptic equations.
The Immersed Interface Method - A Numerical Approach for Partial Differential Equations with Interfaces
Immersed interface methods for one-dimensional moving interface problems
ADI methods for heat equations with discontinuties along an arbitrary interface
Math Modeling
On the rapid evaluation of heat potentials on general regions.
Rapid fourth order accurate methods for the solution of the stokes problem in the presence of an immersed boundary.
The fast solution of Poisson's and the biharmonic equations on irregular regions
Fast parallel iterative solution of Poisson's and the biharmonic equations on irregular regions
An implicit numerical method for fluid dynamics problems with immersed elastic boundaries
Axisymmetric bubbles or drops in a uniform flow
Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations
An adaptive Cartesian grid method for unsteady compressible flow in complex geometries.
Numerical analysis of blood flow in the heart
Lectures on mathematical aspects of physiology
Modeling prosthetic heart valves for numerical analysis of blood flow in the heart
Improved Volume
Boundary Integral and Singularity Methods for Linearized Viscous Flow
On the numerical solution of Helmholtz's equation by the capacitance matrix method
An alternative to unstructured grids for computing gas dynamic flow around arbitrarily complex two-dimensional bodies
TRANAIR: a computer code for transonic analyses of arbitrary configurations.
Dynamics of drop deformation and breakup in viscous flows
A level set approach for computing solutions to incompressible two-phase flow
in Studies in Numerical Analysis
Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods
Computations of multi-fluid flows
A front-tracking method for viscous
A method for solving the transonic full-potential equations for general config- urations
--TR

--CTR
Peter Smereka, The numerical approximation of a delta function with application to level set methods, Journal of Computational Physics, v.211 n.1, p.77-90, 1 January 2006
Anna-Karin Tornberg , Bjrn Engquist, Numerical approximations of singular source terms in differential equations, Journal of Computational Physics, v.200 n.2, p.462-488, November 2004
M. Oevermann , R. Klein, A Cartesian grid finite volume method for elliptic equations with variable coefficients and embedded interfaces, Journal of Computational Physics, v.219 n.2, p.749-769, December, 2006
Shaozhong Deng , Kazufumi Ito , Zhilin Li, Three-dimensional elliptic solvers for interface problems and applications, Journal of Computational Physics, v.184 n.1, p.215-243,
Bertrand Maury, A Fat Boundary Method for the Poisson Problem in a Domain with Holes, Journal of Scientific Computing, v.16 n.3, p.319-339, September 2001
Boyce E. Griffith , Charles S. Peskin, On the order of accuracy of the immersed boundary method: higher order convergence rates for sufficiently smooth problems, Journal of Computational Physics, v.208 n.1, p.75-105, 1 September 2005
Tao Lin , Yanping Lin , Robert Rogers , M. Lynne Ryan, A rectangular immersed finite element space for interface problems, Scientific computing and applications, Nova Science Publishers, Inc., Commack, NY, 2001
Tetsuro Yamamoto, Convergence of consistent and inconsistent finite difference schemes and an acceleration technique, Journal of Computational and Applied Mathematics, v.140 n.1-2, p.849-866, 1 March 2002
Ming-Chih Lai , Zhilin Li , Xiaobiao Lin, Fast solvers for 3D Poisson equations involving interfaces in a finite or the infinite domain, Journal of Computational and Applied Mathematics, v.191 n.1, p.106-125, 15 June 2006
Shi Jin , Xuelei Wang, Robust numerical simulation of porosity evolution in chemical vapor infiltration: II. Two-dimensional anisotropic fronts, Journal of Computational Physics, v.179 n.2, p.557-577, July 2002
Changhoon Lee, Stability characteristics of the virtual boundary method in three-dimensional applications, Journal of Computational Physics, v.184 n.2, p.559-591, 20 January
Elijah P. Newren , Aaron L. Fogelson , Robert D. Guy , Robert M. Kirby, Unconditionally stable discretizations of the immersed boundary equations, Journal of Computational Physics, v.222 n.2, p.702-719, March, 2007
Bruno Lombard , Jol Piraux, Numerical treatment of two-dimensional interfaces for acoustic and elastic waves, Journal of Computational Physics, v.195 n.1, p.90-116, 20 March 2004
Shuwang Li , John S. Lowengrub , Perry H. Leo, A rescaling scheme with application to the long-time simulation of viscous fingering in a Hele-Shaw cell, Journal of Computational Physics, v.225 n.1, p.554-567, July, 2007
Sheng Xu , Z. Jane Wang, An immersed interface method for simulating the interaction of a fluid with moving boundaries, Journal of Computational Physics, v.216 August 2006
D. V. Le , B. C. Khoo , J. Peraire, An immersed interface method for viscous incompressible flows involving rigid and flexible boundaries, Journal of Computational Physics, v.220 n.1, p.109-138, December, 2006
X. Sheldon Wang, An iterative matrix-free method in implicit immersed boundary/continuum methods, Computers and Structures, v.85 n.11-14, p.739-748, June, 2007
Xiaolin Zhong, A new high-order immersed interface method for solving elliptic equations with imbedded interface of discontinuity, Journal of Computational Physics, v.225 n.1, p.1066-1099, July, 2007
Tsung-Min Hwang , Wen-Wei Lin , Wei-Cheng Wang , Weichung Wang, Numerical simulation of three dimensional pyramid quantum dot, Journal of Computational Physics, v.196 n.1, p.208-232, 1 May 2004
Ashley J. James , John Lowengrub, A surfactant-conserving volume-of-fluid method for interfacial flows with insoluble surfactant, Journal of Computational Physics, v.201 December 2004
Donna Calhoun, A Cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions, Journal of Computational Physics, v.176 n.2, p.231-275, March 1, 2002
Anvar Gilmanov , Fotis Sotiropoulos, A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies, Journal of Computational Physics, v.207 August 2005
Jian-Jun Xu , Zhilin Li , John Lowengrub , Hongkai Zhao, A level-set method for interfacial flows with surfactant, Journal of Computational Physics, v.212 n.2, p.590-616, 1 March 2006
George Biros , Lexing Ying , Denis Zorin, A fast solver for the Stokes equations with distributed forces in complex geometries, Journal of Computational Physics, v.193 n.1, p.317-348, January 2004
Guillaume Perigaud , Richard Saurel, A compressible flow model with capillary effects, Journal of Computational Physics, v.209 n.1, p.139-178, 10 October 2005
Boyce E. Griffith , Richard D. Hornung , David M. McQueen , Charles S. Peskin, An adaptive, formally second order accurate version of the immersed boundary method, Journal of Computational Physics, v.223 n.1, p.10-49, April, 2007
bubbles;stokes flow;discontinuous coefficients;cartesian grids;creeping flow;interface tracking;immersed interface methods