A Cartesian Grid Projection Method for the Incompressible Euler Equations in Complex Geometries.
Many problems in fluid dynamics require the representation of complicated internal or external boundaries of the flow.  Here we present a method for calculating time-dependent incompressible inviscid flow which combines a projection method with a "Cartesian grid" approach for representing geometry.  In this approach, the body is represented as an interface embedded in a regular Cartesian mesh. The advection step is based on a Cartesian grid algorithm for compressible flow, in which the discretization of the body near the flow uses a volume-of-fluid representation. A redistribution procedure is used to eliminate time-step restrictions due to small cells where the boundary intersects the mesh. The projection step uses an approximate projection based on a Cartesian grid method for potential flow. The method incorporates knowledge of the body through volume and area fractions along with certain other integrals over the mixed cells. Convergence results are given for the projection itself and for the time-dependent algorithm in two dimensions. The method is also demonstrated on flow past a half-cylinder with vortex shedding.
Introduction
In this paper, we present a numerical method for solving the unsteady incompressible Euler
equations in domains with irregular boundaries. The underlying discretization method is
a projection method [22, 5]. Discretizations of the nonlinear convective terms and lagged
pressure gradient are first used to construct an approximate update to the velocity field;
the divergence constraint is subsequently imposed to define the velocity and pressure at
the new time. The irregular boundary is represented using the Cartesian mesh approach
[58], i.e. by intersecting the boundary with a uniform Cartesian grid, with irregular cells
appearing only adjacent to the boundary. The extension of the basic projection methodology
to the Cartesian grid setting exploits the separation of hyperbolic and elliptic terms of the
method in [5] to allow us to use previous work on discretization of hyperbolic and elliptic
PDE's on Cartesian grids. The treatment of the hyperbolic terms is based on algorithms
developed for gas dynamics, and closely follows the algorithm of Pember et al. [56, 57]. The
Cartesian grid projection uses the techniques developed by Young et al [72] for full potential
transonic flow to discretize the elliptic equation that is used to enforce the incompressibility
constraint. The overall design goals of the method are to be able to use the high-resolution
finite difference methods based on higher-order Godunov methods for the advective terms
in the presence of irregular boundaries that effectively add a potential flow component to
the solution.
Cartesian grid methods were first used by Purvis and Burkhalter [58] for solving the
equations of transonic potential flow; see also [71, 40, 62, 72]. Clarke et al. [23] extended
the methodology to steady compressible flow; see also [30, 29, 28, 52]. Zeeuw and Powell
[74], Coirier and Powell [24], Coirier [25], Melton et al. [51], and Aftosmis et al. [1] have
developed adaptive methods for the steady Euler and Navier-Stokes equations.
For time-dependent hyperbolic problems, the primary difficulty in using the Cartesian
grid approach lies in the treatment of the cells created by the intersection of the irregular
boundary with the uniform mesh. There are no restrictions on how the boundary intersects
the Cartesian grid (unlike the "stair step" approach which defines the body as aligned
with cell edges), and as a result cells with arbitrarily small volumes can be created. A
standard finite-volume approach using conservative differencing requires division by the
volume of each cell; this is unstable unless the time step is reduced proportionally to the
volume. Although in the projection method convective differencing is used for the hyperbolic
terms in the momentum equation, we base our methodology for incompressible flow on the
experience gained for compressible flow in the handling of small cells. (In addition, we will
wish to update other quantities conservatively.) The major issues, then, in designing such
a method are how to maintain accuracy, stability, and conservation in the irregular cells at
the fluid-body interface while using a time step restricted by CFL considerations using the
regular cell size alone. We refer to this as the "small cell problem."
Noh [54] did early work in this area in which he used both cell merging techniques
and redistribution. LeVeque [42, 43] and Berger and LeVeque [12] have developed explicit
methods which use the large time step approach developed by LeVeque [41] to overcome
the small cell problem. Berger and LeVeque [13, 14] have also studied approaches in which
the small cell problem is avoided by the use of a rotated difference scheme in the cells cut
by the fluid-body interface. Both the large time step and the rotated difference schemes
are globally second-order and better than first-order but less than second-order at the fluid-
body interface. Cell merging techniques were also used by Chiang et al. [20], Bayyuk et al.
[4], and Quirk [60, 59] to overcome the small cell problem. Results for this method suggest
it is globally second order accurate and first order at the boundary. Two other approaches
for unsteady compressible flow are based on flux-vector splitting (Choi and Grossman [21])
and state-vector splitting ( -
Oks-uzo-glu [55], Gooch and -
Oks-uzo-glu [31]), but neither of these
approaches avoids the small cell problem.
The advection step of the method presented here uses a different approach for handling
irregular cells based on ideas previously developed for shock tracking by Chern and Colella
[17] and Bell, Colella, and Welcome [8], and extended by Pember et al. [56, 57]. In this approach
the boundary is viewed as a "tracked front" in a regular Cartesian grid with the fluid
dynamics at the boundary governed by the boundary conditions of a stationary reflecting
wall. The basic integration scheme for hyperbolic problems consists of two steps. In the
first step, a reference state is computed using fluxes generated by a higher-order Godunov
method in which the fluid-body boundary is essentially "ignored". In the second step, a
correction is computed to the state in each irregular cell. A stable, but nonconservative,
portion of this correction is applied to the irregular cell. Conservation is then maintained by
a variation of the algebraic redistribution algorithm in [17] which distributes the remainder
of the correction to those regular and irregular cells that are immediate neighbors of the
irregular cell. This redistribution procedure allows the scheme to use time steps computed
from CFL considerations on the uniform grid alone. We adapt this scheme for the advection
step of the projection method.
The projection step requires solution of an elliptic equation on an irregular grid. The
finite-element-based projection developed by Almgren et al. [3] for a regular grid is extended
here to accommodate embedded boundaries using the techniques developed by Young et
al. [72]. The approach of Young et al. defines the approximating space using standard
bilinear elements on an enlarged domain that includes nodes of mixed cells lying outside
the domain. Quadratures in the weak form are restricted to the actual domain to define
the discretization. The resulting linear system can be readily solved with a straightforward
extension of multigrid.
A second-order Cartesian grid projection method has been developed recently by Tau
[67] for the incompressible Navier-Stokes equations. In this formulation velocities are defined
on a staggered grid (which is incompatible with the cell-centered form of the Godunov
method), and the tangential velocity must vanish at the boundaries (i.e. the formulation
is not extensible to the Euler equations). Comments and results in [67] indicate that the
accuracy of the method again is in general first-order at the boundary but globally second-
order. No mention is made of the small cell problem for advecting quantities other than
momentum.
In addition to the approach taken here, there are two other basic approaches to the
treatment of irregular domains: and locally body-fitted grids obtained from mesh map-
pings, and unstructured grids. For body-fitted structured or block-structured grids for
compressible flow, the literature is quite extensive. Typical examples include [37, 35, 16,
11, 27, 7, 65, 68, 9, 61, 73]. The primary objection to mapped grids is the difficulty in
generating them for complex geometries, particularly in three dimensions. There have been
a number of approaches to making mapped grids more flexible, such as using combinations
of of unstructured and structured grids [69, 70, 49, 64, 53, 36, 48], multiblock-structured
grids [2, 66, 63] and overlapping composite grids [10, 19, 34]. On the topic of unstructured
grids, we refer the reader to [39, 38, 45, 46, 50] for compressible flow. L-ohner et al [44] have
developed an unstructured finite element algorithm for high-Reynolds number and inviscid
incompressible flows.
In comparison with the other approaches discussed above, the advantages of the Cartesian
grid approach include not requiring special grid generation techniques to fit arbitrary
boundaries, the ability to use a standard time-stepping algorithm at all cells with no additional
work other than at cells at or near the boundary, and the ability to incorporate
efficient solvers for the elliptic equation. The most serious disadvantage is a loss of accuracy
at the boundary (numerical results show a reduction to first-order accuracy at the
boundary). The loss of second-order accuracy at the boundary indicates that the Cartesian
grid approach is not ideal for calculations designed to resolve viscous boundary layers, but
should be satisfactory for flows where the primary effect of the boundary is to introduce a
potential flow component to the velocity corresponding to the effect of the no-normal flow
condition at the boundary. Problems where this is the case include flows in large containers,
such as utility burners and boilers, and atmospheric flows over irregular orography.
In the next section we review the basic fractional step algorithm, and introduce the
notation of the Cartesian grid method. The subsequent two sections contain descriptions of
the advection step and the projection step, respectively, for flows with embedded boundaries.
In the final two sections we present numerical results and conclusions. All results and
detailed discussion will be for two spatial dimensions; the extension to three dimensions is
straightforward and will be presented in later work.
2 Basic Algorithm
2.1 Overview of Fractional Step Formulation
The incompressible Euler equations written in convective form are
and
r
Alternatively, (2.1) could be written in conservative form as
The projection method is a fractional step scheme for solving these equations, composed
of an advection step followed by a projection. In the advection step for cells entirely in the
flow domain we solve the discretization of (2.1),
\Deltat
for the intermediate velocity U   : For small cells adjoining the body we modify the velocity
update using the conservative formulation of the nonlinear terms (see section 3). The
pressure gradient at t n\Gamma 1
was computed in the previous time step and is treated as a source
term in (2.4). The advection terms in (2.4), namely [(U \Delta r)U
2 , are approximated at
time
2 to second-order in space and time using an explicit predictor-corrector scheme;
their construction is described in section 3.
The velocity field U   is not, in general, divergence-free. The projection step of the
algorithm decomposes the result of the first step into a discrete gradient of a scalar potential
and an approximately divergence-free vector field. They correspond to the update to the
pressure gradient and the update to the velocity, respectively. In particular, if P represents
the projection operator then
\Deltat
\Deltat
\Deltat
Note that the pressure gradient is defined at the same time as the time derivative of velocity,
and therefore at half-time levels.
2.2 Notation
We first introduce the notation used to describe how the body intersects the computational
domain. The volume fraction   i;j for each cell is defined as the fraction of the computational
cell B i;j that is inside the flow domain. The area fractions a i+ 1
2 ;j and a i;j+ 1specify the
fractions of the (i
respectively, that lie inside the flow domain.
We label a cell entirely within the fluid ( as a full cell or fluid cell, a cell entirely
outside of the flow domain ( as a body cell, and a cell partially in the fluid (0 !
as a mixed cell. A small cell is a mixed cell with a small volume fraction. We note that it
is possible for the area fraction of an edge of a fluid cell to be less than one if the fluid cell
abuts a body or mixed cell.
In the method presented here, the state of the fluid at time t n is defined by U n
ij ); the velocity field in cell B i;j at time t n ; and p n\Gamma
the pressure at node (i\Gamma
2 . Gp n\Gamma 1ij is the pressure gradient in cell B i;j at time t
For the construction
of the nonlinear advective terms at t n+ 1
velocities are also defined on all edges of full and
mixed cells at t n+ 1
process requires values of the velocity and pressure gradients in the
cells on either side of an edge at t n : We must therefore define, at each time step, extended
states in the body cells adjoining mixed or full cells. We do this in a volume-weighted
fashion:
U ext
The pressure gradient is extended in the same manner. Here nbhd(B i;j ), defined as the
neighborhood of B i;j ; refers to the eight cells (in two dimensions) that share an edge or
corner with
3 Advection step
3.1 Momentum
The algorithm is a predictor-corrector method, similar to that used in [5], but with some
modifications as discussed in [6]. The details of the current version without geometry are
given in [3]. For simplicity we will assume that the normal velocity on the embedded
boundary is zero; the treatment of a more general Dirichlet boundary condition such as
inflow is straightforward.
In the predictor we extrapolate the velocity to the cell edges at t n+ 1
using a second-order
Taylor series expansion in space and time. The time derivative is replaced using (2.1).
For edge (i
U n;lim
extrapolated from B i;j , and
U n;lim
extrapolated from
Analogous formulae are used to predict values at each of the other edges of the cell. In
evaluating these terms the first-order derivatives normal to the edge (in this case U n;lim
x are
evaluated using a monotonicity-limited fourth-order slope approximation [26]. The limiting
is done on the components of the velocity individually, with modifications as discussed below
in cells near the body.
The transverse derivative terms ( d
vU y in this case) are evaluated as in [6], by first extrapolating
from above and below to construct edge states, using normal derivatives only,
and then choosing between these states using the upwinding procedure defined below. In
particular, we define
\Deltay
U n;lim
U n;lim
where U n;lim
y are limited slopes in the y direction, with similar formulae for the lower edge
of B i;j . In this upwinding procedure we first define the normal advective velocity on the
edge:
(We suppress the
spatial indices on bottom and top states here and in the next
equation.) We now upwind b
U based on b
U i;j+ 1=
After constructing b
similar manner, we use these upwind values to form the
transverse derivative in (3.1):
2\Deltay
(b v adv
U
We use a similar upwinding procedure to choose the appropriate states U i+ 1
given the left
and right states, U n+ 1;L
U
We follow a similar procedure to construct U i\Gamma 1
In general the normal velocities at the edges are not divergence-free. In order to make
these velocities divergence-free we apply a MAC projection [6] before the construction of
the convective derivatives. The equation
solved for OE in all full or mixed cells, with
a
a i;j+ 1v n+ 1i;j+ 1\Gamma a
\Deltax
\Deltay
Note that in regions of full cells the D MAC G MAC stencil is simply a standard five-point
cell-centered stencil for the Laplacian. For edges within the body that are needed for the
corrector step below, the MAC gradient is linearly extrapolated from edges that lie partially
or fully in the fluid. For example, if edge (i
lies completely within the body, but
edges
lie at least partially within the fluid, then we define
While it is possible this could lead to instabilities, none have yet been observed in numerical
testing. We solve this linear system using standard multigrid methods (see [15]), specifically,
V-cycles with Jacobi relaxation.
In the corrector step, we form an approximation to the convective derivatives in (2.1)
2\Deltay
(v MAC
where
The intermediate velocity U   at time n+ 1 is then defined on all full and mixed cells as
The upwind method is an explicit difference scheme and, as such, requires a time-step
restriction. When all cells are full cells, a linear, constant-coefficient analysis shows that for
stability we require
\Deltay
where oe is the CFL number. The time-step restriction of the upwind method is used to set
the time step for the overall algorithm.
As mentioned earlier, for incompressible flow the two forms of the momentum equation,
and (2.3), are analytically equivalent. For flows without geometry we have successfully
used the convective difference form of the equations, (2.4). By contrast, we could use the
conservative form of the equations,
\Deltat +r \Delta F n+ 1
2 ); and evaluate the conservative update as
R
In the presence of embedded boundaries, the convective update form of the equations is stable
even for very small cells, but this update is calculated as if the cell were a full cell, i.e. it
doesn't ``see'' the body other than through the MAC-projected normal advection velocities
(and the modification of the limited slopes as discussed at the end of this section). To
compute the conservative update, however, one ignores entirely the portion of the cell not
in the fluid, and integrates the flux only along the parts of the edges of the cell that lie
within the fluid. Thus the presence of the body is correctly accounted for, but using the
conservative form requires that the time step approach zero as the cell volume goes to zero,
which is the small-cell time step restriction which we seek to avoid.
A solution in this case is to use a weighted average of the convective and conservative
updates, effectively allowing as much momentum to pass into a small cell in a time step
as will keep the scheme stable. The momentum that does not pass into the small cell is
redistributed to the neighboring cells, to maintain conservation in the advection step. This
approach is modeled on the algorithm of [57, 56], and is based on the algebraic redistribution
scheme of Chern and Colella [18]. The algorithm is as follows:
(1) First, in all cells construct the reference state, e
U   , defined as
e
using the advective algorithm described for full cells.
(2) On mixed cells only construct an alternative e e
U

using a conservative
U

(a
(a i;j+ 1F y
\Deltay
where the fluxes are defined as
This solution enforces no-flow across the boundary of the body, but is not necessarily stable.
Define the difference ffiM
U

U
(3) The conservative solution can be written e e
U

U
; however, this solution is
not stable for   - 1: Define instead b
U
U
ffiM i;j on mixed cells, b
U
U
i;j on full
cells. This allows the mixed cell state to keep the fraction of ffiM which will keep the scheme
stable given that the time step is set by the full-cell CFL constraint.
Now redistribute the remaining fraction of ffiM from each mixed cell, the amount
to the fluid and mixed cells among the eight neighbors of B i;j in a volume-weighted
fashion. Since we redistribute the extensive rather than intensive quantity (e.g.
momentum rather than momentum density), the resulting redistribution has the form
c c
U
U   i;j
where
(5) Subtract the pressure gradient term from the solution for all full and mixed cells, treating
it as a source term:
U
A second modification to the algorithm for cells at or near the boundary is to the slope
calculations; this is based on the principle that no state information from cells entirely
within the body should be used in calculating slopes. Despite the use of extended states in
the creation of edge states, the slope calculation uses state information only from full and
mixed cells. In a mixed or full cell for which the fourth-order stencil would require using
an extended state, the slope calculation reduces to the second-order monotonicity-limited
formula. If the second-order formula would require an extended state, then the slope is set
to zero. This has the effect of adding dissipation to the scheme; if the slopes throughout
the domain were all set to zero the method would reduce to first-order.
Two comments about the algorithm are in order here. First, we note that while the use
of extended states seems to be a low-order approach to representing the body, numerical
results show that it is adequate to maintain first-order accuracy at the boundary. An
alternative to extended states is the so-called "thin wall approximation" of [57]; with this
approximation extended states are not defined, rather upwinding always chooses the values
coming from the fluid side (rather than the body side) of the fluid-body interface. In our
numerical testing, however, we have found that the extended states version of the algorithm
performs slightly better.
Second, we are unable to provide a proof that the redistribution algorithm is stable,
but in extensive numerical testing of the redistribution algorithm for compressible flow (see
[57, 32, 47]) no stability problems have been observed, and none have been observed in our
incompressible flow calculations. It is clear that the algorithm removes the original small-cell
stability problem by replacing division by arbitrarily small volume fractions with division
by one in step (3). The redistribution algorithm is such that the amount redistributed into
a cell is proportional to the volume of that cell; this feature combined with the observation
in our calculations that the values of ffiM are small relative to the values of U in the mixed
cells gives a heuristic argument as to why redistribution would be stable.
3.2 Passive Scalars
The algorithm for advecting momentum extends naturally to linear scalar advection. In
the third numerical example shown in this paper, a passive scalar enters the domain at
the inflow boundary and is advected around the body. We now briefly describe the scalar
advection routine, assuming the scalar s is a conserved quantity, i.e.
First, at each time step extended states are defined for the passive scalar just as for
velocity. In the predictor we extrapolate the scalar to the cell edges at t n+ 1
using a second-order
Taylor series expansion in space and time. For edge (i
\Deltat
extrapolated from B i;j , and
s n+ 1;R
\Deltat
extrapolated from
Analogous formulae are used to predict values at each of the other edges of the cell.
As with velocity, derivatives normal to the edge (in this case s n;lim
x are evaluated using
the monotonicity-limited fourth-order slope approximation [26], with modifications due to
geometry.
The transverse derivative terms (vc s y in this case) are defined using c
where
and
We use a similar procedure to choose s n+ 1i+ 1;j given s n+ 1
if u MAC
s n+ 1;R
if u MAC
In the corrector step for full cells, we update s conservatively:
s
\Deltat
\Deltax
\Deltat
\Deltay (v MAC
In the corrector step for mixed cells, we follow the three steps below:
(1) Construct the reference state e
s   using the convective formulation:
\Deltat
\Deltat
2\Deltay (v MAC
(2) Construct e e s
using the conservative formulation:
(a i+ 1;j u MAC
\Deltax
(a i;j+ 1u MAC
\Deltay
Define the difference ffis
s   i;j
Construct b s
s
Finally, redistribute the remaining fraction of ffis from each mixed cell, the amount
to the full and mixed cells among the eight neighbors of B i;j in a volume-weighted
fashion:
In full cells not adjacent to any mixed cells, set s n+1
Projection Step
Since U   is defined on cell centers, the projection used to enforce incompressibility at time
must include a divergence operator that acts on cell-centered quantities, unlike the
MAC projection. The projection we use here is approximate; i.e. P 2 6= P: The operator, as
well as the motivation for using an approximate rather than exact projection, is described
in detail in [3]. Here we outline the algorithm for the case of no-flow physical boundaries
with no geometry, and for the case of embedded boundaries. The basic approach for the
embedded boundaries utilizes the same discretization as was used by Young et al [72] for
full potential transonic flow. The necessary modifications for inflow-outflow conditions are
described in [3] and their use is demonstrated in the numerical examples.
This projection is based on a finite element formulation. In particular, we consider the
scalar pressure field to be a C 0 function that is bilinear over each cell; i.e., the pressure is
in
s (x) is the space of polynomials of degree t in the x direction on each cell with C s
continuity at x-edges. For the velocity space we define
in x and a discontinuous linear function of y in each cell, with a similar form for v: For
mixed cells, we think of the representations as only being defined on the portion of each cell
within the flow
domain,\Omega ; although the spaces implicitly define an extension of the solution
over the entire cell. Note that in the integral used to define the projection (see (4.3) below)
the domain of integration is limited to the actual flow domain; we do not integrate over the
portion of mixed cells outside the domain.
For use in the predictor and corrector, the velocity and pressure gradient are considered
to be average values over each cell. The vector space, V h , contains additional functions that
represent the linear variation within each cell. These additional degrees of freedom make
large enough to contain rOE for OE 2 S h . We establish a correspondence between these
two representations by introducing an orthogonal decomposition of V h . In particular, for
each we define a piecewise constant component V and the variation V
so that for each cell B i;j ;
R
construction these two components are
orthogonal in L 2 so they can be used to define a decomposition of V h into two subspaces
represent the cell averages and the orthogonal linear variation, respec-
tively. The decomposition of V h induces a decomposition of rOE for all OE 2 S h ; namely,
We now define a weak form of the projection on V h , based on a weak divergence on V h .
In particular, we define a vector field V d in V h to be divergence-free in the
domain\Omega if
Z\Omega
Using the definition (4.2) we can then project any vector field V into a gradient rOE and
weakly divergence-free field V d (with vanishing normal velocities on boundaries) by solving
Z
Z
for
setting V Here we define the /'s to be
the standard basis functions for S h , namely / i+ 1;j+ 1(x) is the piecewise bilinear function
having node values / i+1=2;j+1=2
For the purposes of the fractional step scheme we wish to decompose the vector field
\Deltat
into its approximately divergence free-part
\Deltat
and the update to the pressure,
Since the finite-difference advection scheme is designed to handle cell-based quantities
that are considered to be average values over each cell, the quantity (GOE) ? is discarded at
the end of the projection step. This makes the projection approximate, i.e. if D(V
then in practice, we solve the system defined by (4.3) (again
using standard multigrid techniques with V-cycles and Jacobi relaxation), and set
Z
as the approximation to U n+1 \GammaU n
\Deltat in (2.5). (See [3] for more details.)
The left-hand-side of equation (4.3) is, in discrete form, a nine-point stencil approximating
the Laplacian of OE and the right hand side, for is a standard four-point
divergence stencil. Without embedded boundaries, the method reduces to constant coefficient
difference stencils for divergence, gradient, and the Laplacian operator. These stencils
are, for
(DV
and, letting
(DGOE)
Body cells contribute nothing to the integrals in (4.3); in mixed cells the integrals are
computed only over the portion of each cell that lies in the fluid. This calculation can
be optimized by pre-computing the following integrals in each mixed cell; these integrands
result from the products of the gradients of the basis functions in (4.3). We define
I
Z
I
Z
I
Z
"\Omega
are functions defined on each cell such that x 0 at the center of
\Sigma1=2 at the left and right edges of B i;j ; and y \Sigma1=2 at the top and bottom
edges of B
Linear combinations of these integrals form the coefficients used in the divergence, gradient
and Laplacian operators. For example, the divergence of a barred vector becomes
(DV
I a
I a
I a
I a
)=\Deltay
If we let (GOE) full be (GOE) as it would be defined if the cell were a full cell, then in a
mixed cell the average of the gradient over the fluid portion of the cell is
full
I b
full
I a
where, for
The full Laplacian operator requires the third integral I c because of the quadratic terms
which results from the inner products of OE i+ 1;j+ 1and OE k+1=2;'+1=2 ; we leave the derivation
to the reader. Note that for a full cell  and so the stencils above for
mixed cells reduce to those for full cells.
Note that solving (4.3) defines OE i+ 1;j+ 1at each node (i
which the support
of
intersects the fluid region. Thus, the pressure is defined even at nodes
contained in the body region, as long as they are within one cell of the fluid region.
5 Numerical Results
In this section we present results of calculations done using the Cartesian grid representation
of bodies and/or boundaries of the domain. The first two sets of results are convergence
studies; the first tests the accuracy of the projection alone, the second tests the accuracy of
the full method for the Euler equations. In both cases, the results are presented in the form
of tables which show the norms of the errors, as well as the calculated rates of convergence.
The error for a calculation on a grid of spacing h is defined as the difference between the
solution on that grid and the averaged solution from the same calculation on a grid of
spacing h=2: In the first convergence study, the column 128-256 refers to the errors of the
solution on the 128 2 grid as calculated by comparing the solution on the 128 2 grid to the
solution on the 256 2 grid; similarly for 256-512. Once the errors are computed pointwise for
each calculation, the L 1 , L 2 , and L1 norms of the errors are calculated. The convergence
rates of the method can be estimated by taking the log 2 of the ratio of these norms. This
provides a heuristic estimate of the convergence rate assuming that the method is operating
in its asymptotic range. We present two separate measures of the error. In the first we
compute the norms over the entire domain ("All cells"). In the second we examine the
error in an interior subdomain in order to measure errors away from the boundary. For
the subdomain we have selected the region covered by full cells in the 128 2 grid ("Full 128
cells").
The first convergence study tests the accuracy of the projection alone. The domain is
the unit square, and three bodies are placed in the interior. A circle is centered at (:75; :75)
and has radius .1; ellipses are centered at (:25; :625) and (:5; :25) and have axes (:15; :1) and
respectively. The boundary conditions are inflow at
and no-flow boundaries at
In this study, the initial data inside the numerical domain is (u; this data is
then projected to define a velocity field which is approximately divergence-free in the region
of the domain not covered by the bodies. The pressure gradient which is used to correct
the initial data represents the deviation of the potential flow solution from uniform flow.
In

Tables

1 and 2 we present the results of this convergence study. The almost-second-
order convergence in the L 1 norm, first-order convergence in the L1 norm, and approximately
h 1:5 convergence in the L 2 norm, are what we would expect of a solution which is
second-order in most of the domain and first-order at boundaries. This is consistent with
the observation that the maximum absolute error on the mixed cells is an order of magnitude
higher than the maximum error on the full cells for the finer grids. Figure 1 shows a
contour plot of the magnitude of the error of the 256 2 calculation.
The second convergence study is of flow through a diverging channel. In this case we
evaluate the velocity field at time in order to demonstrate the order of the complete
algorithm for flow that is smooth in and near the mixed cells. The problem domain is 4 x
1, and the fluid is restricted to flow between the curves y top and y bot ; defined as
3:
3: - x - 4:
and
with The calculation is run at CFL number 0.9, and the flow is
initialized with the potential flow field corresponding to inflow velocity 1: at the left
edge, as computed by the initial projection. Tables 3 and 4 show the errors and convergence
rates for the velocity components; here 128-256 refers to the errors of the solution on the
128x32 grid as calculated by comparing the solution on the 128x32 grid to the solution
on the 256x64 grid; similarly for 256-512. As before, we compute the errors both on the
full domain and on a subdomain defined as the region covered by full cells on the 128x32
grid. (In this case, the norms are scaled, as appropriate, by the area of the domain.) Again
we see rates corresponding to global second-order accuracy but first-order accuracy near
the boundaries. Figures 2a-b are contour plots of the log 10 of the error in the velocity
components. The error is clearly concentrated along the fluid-body boundaries; for the sake
of the figures we have defined all errors less than 1% of the maximum error in each (including
the values in the body cells) to be equal to 1% of the maximum error. There are ten contour
intervals in each figure spanning these two orders of magnitude. In time the error advects
further along but does not contaminate the interior flow. We note
that for this particular flow the velocity field near the boundary is essentially parallel to the
boundary, and in a more general case we might expect to see more contamination of the
interior flow. Note, however, that the maximum error is less than :3% of the magnitude of
the velocity.
The third example we present is that of flow past a half-cylinder. There is an extensive
experimental and computational literature on the subject of flow past a cylinder in an
infinite domain at low to moderate Reynolds number; see, for example [33] and [61] for recent
experimental and computational results, respectively. Since the methodology presented here
is for inviscid flow, we compute flow past a half-cylinder rather than a full cylinder so as to
force the separation point to occur at the trailing edge. However, we present this calculation
to demonstrate the type of application for which the Cartesian grid methodology would be
most useful: a calculation in which one might be interested most in the flow features away
from the body (i.e. the shed vortices downstream from the cylinder), but which requires
the presence of the body in order to generate or modify those features.
The resolution of our calculation is 256 x 64, the domain is 4 x 1, and the diameter of the
half-cylinder is .25. The inflow velocity at the left edge is the boundary conditions are
outflow at the right edge, no-flow boundaries at the top and bottom. The initial conditions
are the potential flow with uniform inflow as calculated from the initial projection, combined
with a small vortical perturbation to break the symmetry of the problem. Shown here are
"snapshots" of the vorticity (Figure 3a) and a passively advected scalar (Figure 3b) at late
enough time that the flow is periodic, and that the perturbation has been advected through
the domain. The scalar was advected in from the center of the inflow edge. The Strouhal
number is calculated to be approximately D=(U1T
is the observed period of vortex shedding, is the cylinder diameter, and
is the free-stream (i.e. inflow) velocity. One would expect a value between .2 and .4 (see
[33] and the references cited there), so this seems a reasonable value given the limitations
of the comparison.
All cells Full 128 cells

Table

1: Errors and convergence rates for the x-component of the pressure gradient.
All cells Full 128 cells

Table

2: Errors and convergence rates for the y-component of the pressure gradient.
All cells Full 128 cells
L1 2.36e-3 .80 1.35e-3 1.28e-3 1.08 6.04e-4

Table

3: Errors and convergence rates for the x-component of the velocity.
All cells Full 128 cells

Table

4: Errors and convergence rates for the y-component of the velocity.
6 Conclusion
We have presented a method for calculation of time-dependent incompressible inviscid flow
in a domain with embedded boundaries. This approach combines the basic projection
method, using an approximate projection, with the Cartesian grid representation of ge-
ometry. In this approach, the body is represented as an interface embedded in a regular
Cartesian mesh. The adaptation of the higher-order upwind method to include geometry
is modeled on the Cartesian grid method for compressible flow. The discretization of the
body near the flow uses a volume-of-fluid representation with a redistribution procedure.
The approximate projection incorporates knowledge of the body through volume and area
fractions, and certain integrals over the mixed cells. Convergence results indicate that the
method is first-order at and near the body, and globally second-order away from the body.
The method is demonstrated on flow past a half-cylinder with vortex shedding, and is shown
to give a reasonable Strouhal number.
The method here is presented in two dimensions; the extension to r \Gamma z and three
dimensions and to variable density flows, and the inclusion of this representation with
an adaptive mesh refinement algorithm for incompressible flow are being developed. The
techniques developed here are also being modified for use in more general low Mach number
models. In particular, we are using this methodology to model low Mach number combustion
in realistic industrial burner geometries and to represent terrain in an anelastic model of
the atmosphere.

Acknowledgements

We would like to thank Rick Pember for the many conversations on the subtleties of the
Cartesian grid method.



--R

Adaptation and surface modeling for cartesian mesh methods
Techniques in multiblock domain decomposition and surface grid generation
A numerical method for the incompressible Navier-Stokes equations based on an approximate projection
A simulation technique for 2-d unsteady inviscid flows around arbitrarily moving and deforming bodies of arbitrary geometry
A second-order projection method for the incompressible Navier-Stokes equations
An efficient second-order projection method for viscous incompressible flow
Adaptive mesh refinement on moving quadrilateral grids
Conservative front-tracking for inviscid compressible flow
A projection method for viscous incompressible flow on quadrilateral grids
A flexible grid embedding technique with application to the Euler equations
Automatic adaptive refinement for the Euler equations
An adaptive Cartesian mesh algorithm for the Euler equations in arbit rary geometries



Progress in finite-volume calculations for wind- fuselage combinations
A conservative front tracking method for hyperbolic conservation laws
A conservative front tracking method for hyperbolic conservations laws
Composite overlapping meshes for the solution of partial differential equations
Simulation of unsteady inviscid flow on an adaptively refined Cartesian grid

Numerical solution of the navier-stokes equations
Euler calculations for multielement airfoils using Cartesian grids
An accuracy assessment of Cartesian-mesh approaches for the Euler equ ations

A direct Eulerian MUSCL scheme for gas dynamics

Cartesian Euler methods for arbitrary aircraft configurations
An adaptive multigrid applied to supersonic blunt body flow
Euler calculations for wings using Cartesian grids.

An adaptive multifluid interface-capturing method for compressible flow in complex geometries
An experimental study of the parallel and oblique vortex shedding from circular cylinders
A fourth-order accurate method for the incompressible Navier-Stokes equations on overlapping grids

Inviscid and viscous solutions for airfoil/cascade flows using a locally implicit algorithm on adaptive meshes
Transonic potential flow calculations using conservative form
Improvements to the aircraft Euler method.
Calculation of inviscid transonic flow over a complete aircraft.

A large time step generalization of Godunov's method for systems of conservation laws





Multidimensional numerical simulation of a pulse combustor
A hybrid structured-unstructured grid method for unsteady turbomachinery flow compuations
A solution-adaptive hybrid-grid method for the unsteady analysis of turbomachinery
Accurate multigrid solution of the Euler equations on unstructured and adaptive meshes
3d applications of a cartesian grid euler method
A finite difference solution of the Euler equations on non-body fitted grids

Cel: A time-dependent



Prediction of critical Mach number for store config- urations
A Cartesian grid approach with hierarchical refinement for compressible flows

A fractional step solution method for the unsteady incompressible Navier-Stokes equations in generalized coordinate systems


A zonal implicit procedure for hybrid structured-unstructured grids
Computations of unsteady viscous compressible flows using adaptive mesh refinement in curvilinear body-fitted grid systems
A general decomposition algorithm applied to multi-element airfoil grids
A second-order projection method for the incompressible Navier-Stokes equations in arbitrary domains
Boundary fitted coordinate systems for numerical solution of partial differential equations - A review
A three-dimensional struc- tured/unstructured hybrid Navier-Stokes method for turbine blade rows
Mixed structured-unstructured meshes for aerodynamic flow sim- ulation
A method for solving the transonic full-potential equations for general configurations
A locally refined rectangular grid finite element method: Application to computational fluid dynamics and computational physics
fractional step method for time-dependent incompressible Navier-Stokes equations in curvilinear coor- dinates
An adaptively refined Cartesian mesh solver for the Euler equations
--TR

--CTR
Caroline Gatti-Bono , Phillip Colella, An anelastic allspeed projection method for gravitationally stratified flows, Journal of Computational Physics, v.216 August 2006
Yu-Heng Tseng , Joel H. Ferziger, A ghost-cell immersed boundary method for flow in complex geometry, Journal of Computational Physics, v.192 December
Jiun-Der Yu , Shinri Sakai , James Sethian, A coupled quadrilateral grid level set projection method applied to ink jet simulation, Journal of Computational Physics, v.206 n.1, p.227-251, 10 June 2005
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
R. Ghias , R. Mittal , H. Dong, A sharp interface immersed boundary method for compressible viscous flows, Journal of Computational Physics, v.225 n.1, p.528-553, July, 2007
Stphane Popinet, Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries, Journal of Computational Physics, v.190 n.2, p.572-600, 20 September
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
M. P. Kirkpatrick , S. W. Armfield , J. H. Kent, A representation of curved boundaries for the solution of the Navier-Stokes equations on a staggered three-dimensional Cartesian grid, Journal of Computational Physics, v.184 n.1, p.1-36,
S. Marella , S. Krishnan , H. Liu , H. S. Udaykumar, Sharp interface Cartesian grid method I: an easily implemented technique for 3D moving boundary computations, Journal of Computational Physics, v.210 n.1, p.1-31, 20 November 2005
cartesian grid;incompressible Euler equations;projection method