18

me
Congrs Franais de Mcanique Grenoble, 27-31 aot 2007
The LS-STAG Method : A new Immersed Boundary / Level-Set Method for
the Computation of Incompressible Viscous Flows in Complex Geometries
Yoann Cheny & Olivier Botella
Nancy Universits
LEMTA - UMR 7563 (CNRS-INPL-UHP)
2, avenue de la Fort de Haye, B.P. 160, 54504 Vanduvre-ls-Nancy (France)
[email protected] ; [email protected]
Abstract :
This paper concerns the development of a new Cartesian grid / immersed boundary (IB) method for the computa-
tion of incompressible viscous ows in irregular geometries. In IB methods, the irregular boundary is not aligned
with the computational grid, and of upmost importance for accuracy and stability is the discretization in cells
which are cut by the boundary, the so-called cut-cells. In this paper, we present a progress report on a new
nite volume discretization where the irregular boundary is implicitly represented by its signed distance function
(the level-set function). With the help of level-set calculus tools, we are able to preserve the Cartesian structure of
the discretized Navier-Stokes equations in the cut-cells. The resulting discrete systems are efciently solved with
a black box multigrid solver for structured grids. The method is validated for the circular cylinder ow at low
Reynolds number.
Key-words :
Incompressible ows ; Computational uid dynamics ; Immersed boundary methods
1 Introduction
Much attention has recently been devoted to the extension of Cartesian grid ow solvers to
complex geometries by immersed boundary (IB) methods (see [7] for a recent review). In these
methods, the irregular boundary is not aligned with the computational grid, and the treatment of
the cells which are cut by the boundary remains an important issue. Indeed, the discretization
in these cut-cells should be designed such that : (a) the overall accuracy of the method is
not severely diminished and (b) the high computational efciency of the structured solver is
preserved.
Two major classes of IB methods can be distinguished on the basis of their treatment of
cut-cells. Classical IB methods use a nite volume/difference structured solver in Cartesian
cells away from the irregular boundary, and discard the discretization of ow equations in the
cut-cells [7]. Instead, special interpolations are used for setting the value of the dependent
variables in the latter cells. Thus, strict conservation of quantities such as mass and momentum
is not observed near the irregular boundary. Numerous revisions of these interpolations are still
proposed for improving the accuracy and consistency of this class of IB methods [5]. A second
class of IB methods (also called cut-cell methods, see e.g. [3]) aims for actually discretizing
the ow equations in cut-cells. However, the calculation of uxes in cut-cells relies usually on
unstructured techniques, and their negative impact on the computational efciency of the code
is difcult to evaluate.
The purpose of this communication is to present a newIBmethod for incompressible viscous
ows which takes the best aspect of both approaches. As the cut-cell method of Ref. [3], our
method is based on the symmetry preserving nite-volume discretization on Cartesian grids
by Verstappen & Veldman [10]. However, the major difference is that we have undertaken
1
18
me
Congrs Franais de Mcanique Grenoble, 27-31 aot 2007
representing the irregular boundary by its level-set function [8]. With the help of level-set
calculus tools, the uxes in Cartesian and cut-cells are discretized in a unied fashion, ensuring
that the Cartesian structure of the stencils is preserved. Our method has the following distinctive
features :
In contrast to standard IB methods, ow variables are actually computed near the irregular
boundary, and not interpolated [5]. This feature is of upmost importance for computing
the boundary layer of viscous ows. It is even more important for the viscoelastic ow
computations we intend to perform in the future, for which the basic solver we present in
this paper is the building block.
The Cartesian structure of the discrete systems enables the use of efcient black box
solvers for structured grids, resulting in a highly computationally efcient method.
In the nal part of the paper, our method is validated for the circular cylinder ow in the
laminar regime, which has become the standard benchmark test for IB methods.
2 Numerical method
Let be a rectangular computational domain and its surface. The governing equations are
the incompressible Navier-Stokes equations in integral form. In the following, we will consider
the nite-volume discretization of the continuity equation :
_

v n dS = 0, (1)
where v = (u, v) is the velocity, and of the u-momentum equation :
d
dt
_

u dV +F
c
+
_

p e
x
n dS
1
Re
F
d
= 0, (2)
where p is the pressure, Re is the Reynolds number, F
c

(v n) u dS and F
d

u
n dS are the convective and diffusive ux respectively.
Basic discretization on the MAC mesh
The Cartesian method on which our IB method is based is the second-order nite volume dis-
cretization of Verstappen & Veldman [10], which has the ability to preserve on non-uniform
cell distributions the conservation properties (for total mass, momentum and kinetic energy)
of the original MAC method [4]. The staggered arrangement of the discrete velocities on the
MAC mesh is illustrated in Fig. 1. The computational domain is divided in non-uniform Carte-
sian cells
i,j
=

x
i1
, x
i
_

y
j1
, y
j
_
, of size V
ij
= x
i
y
j
, and center x
c
ij
= (x
c
i
, y
c
j
).
The surface of cell
i,j
is divided in 4 elementary faces as
i,j
=
e

w

n

s
using
the usual compass notations. Cell
i,j
is used as the control volume for discretizing the con-
tinuity equation (1), while the staggered control volume
u
i,j
=

x
c
i
, x
c
i+1
_

y
j1
, y
j
_
, of
size V
u
ij
=

x
i
y
j
with

x
i
=
1
2
x
i
+
1
2
x
i+1
, is used for the u-momentum equation (2).
Compass notations are used for intermediate steps of the discretization of (2). For example, we
denote F
c
e
the convective ux through the east face
u
e
of the staggered control volume, i.e. :
F
c
e

_
y
j
y
j1
(v n
e
) u(x
c
i+1
, y) dy. (3)
2
18
me
Congrs Franais de Mcanique Grenoble, 27-31 aot 2007
Figure 1: Location of the variables near the cell
i,j
on the MAC mesh.
The LS-MAC mesh for discretization in irregular domains
We consider now an irregular uid domain
f
which is embedded in the computational domain
. To keep track of the irregular boundary
f
, we employ a signed distance function (x) (i.e. ,
the level set function [8]) such that is negative in the uid region
f
, is positive in the solid
region
s
= \
f
, and such that the boundary
f
corresponds to the zero level-set of this
function, i.e. :
(x)
_
_
_
, x
f
,
0, x
f
,
+, x
s
,
(4)
where is the distance from x to the closest point on the irregular boundary. This leads to the
modication of the MAC mesh that is described in Fig. 2, and that will be subsequently referred
to as the LS-MAC mesh. In this gure, the cell
i,j
which is part uid / part solid, is commonly
called a cut cell. The discretization of the uxes in cut cells is performed with the help of an
additional variable, denoted
i,j
, that stores the values of the level-set function at the upper
right corner of the cells. The level-set values are used to efciently compute quantities relevant
to the irregular boundary and the cut-cells, such as outward normal vectors or surface integrals
over
f
, areas and volumes of the cut-cells, etc. . . In this respect, the most important quantity
for the LS-MAC discretization is the uid portion of the cell faces. For example in Fig. 2, by
using one-dimensional linear interpolation of (x
i
, y) in [y
j1
, y
j
], the length of the face
e
which belongs to the uid domain is :
y
n
y
j1
=
u
i,j
y
j
, with
u
i,j
=

i,j1

i,j1

i,j
since (x
i
, y
n
) = 0.
The quantity
u
i,j

_
0, 1

, which is called the cell-face ratio, shall be used to :

u
i1,j
+
u
i,j
_
x
i
y
j
. (5)
3
18
me
Congrs Franais de Mcanique Grenoble, 27-31 aot 2007
Figure 2: Location of the variables near the cut cell
i,j
on the LS-MAC Mesh. The intersection of the
east face
u
with the irregular boundary is (x
i
, y
n
). See Fig. 1 for complementary notations.
Implement boundary conditions in the fashion of the ghost uid method [8]. For example,
the diffusive ux through the solid north boundary
u
n
in Fig. 2 is computed as :
F
d
n

_

u
n
u
y
dx

=

x
i
u(x
i
, y
n
) u
i,j
1
2

u
i,j
y
j
, (6)
with u(x
i
, y
n
) given by the Dirichlet conditions at the irregular boundary.
These points will be discussed in greater lengths in a forthcoming paper (see also [1]). In
the following, we will mainly detail the discretization of the continuity equation (1).
Discretization of the continuity equation
As in the Cartesian method of Verstappen & Veldman [10], the starting point of the LS-MAC
discretization concerns the continuity equation (1), that reads in cell
i,j
:
u
i,j
u
i1,j
+v
i,j
v
i,j1
= 0, (7)
where mass uxes across cell faces are denoted with a bar ; for example the mass ux through
the east face
u
e
in Fig. 2 is :
u
i,j

_
y
n
y
j1
u(x
i
, y) dy. (8)
In order to easily discretize this integral, we rst locate the discrete unknown u
i,j
in the mid-
dle of the uid part of the face as u
i,j
u(x
i
, y
j1
+
1
2

u
i,j
y
j
). Then, by using midpoint
quadrature, we nally discretize (8) as :
u
i,j

=
u
i,j
y
j
u
i,j
. (9)
After similar calculations for the other cell faces, the discretization of (7) is :
y
j
_

u
i,j
u
i,j

u
i1,j
u
i1,j
_
+x
i
_

v
i,j
v
i,j

v
i,j1
v
i,j1
_
= 0, (10)
and, in the case of a regular Cartesian cell ( = 0 or 1 only), the discrete continuity equation
reduces to the one of the standard MAC discretization.
4
18
me
Congrs Franais de Mcanique Grenoble, 27-31 aot 2007
Further computational details
Discretization of the umomentum equation (2) is performed in staggered control volumes as
shown in Fig. 2. Volume integrals are discretized as :
d
dt
_

u
u dV

= V
u
i,j
du
i,j
dt
, (11)
and the volume of the cut-cells is efciently computed with the cell face ratios (see Eq. (5)).
Discretization of the convective ux follows the lines of the symmetry preserving discretization
of [10]. For example, Eq. (3) is discretized as F
c
e

=
1
2
_
u
i,j
+u
i+1,j
_
u
e
, and central interpolation
with constant weights for u
e
is used. The discretization of the diffusive uxes is more involved
and is not detailed in this paper. The pressure gradient in Eq. (2) is simply computed as :
1
V
u
i,j
_ _

u
e
p dx
_

u
w
p dx
_

=
p
i+1,j
p
i,j
x
p
i
, (12)
with x
p
i
=
1
2
x
i
+
1
2
x
i+1
for the cut-cell of Fig. 2.
The time stepping method we use is the semi-implicit AB/ BDF 2 projection scheme (see
e.g. [2]). With our LS-MAC discretization, the pressure Poisson equation of the projection step
is a symmetric linear system, as in the Cartesian case. This linear system is efciently solved
with a black-box multigrid solver for Cartesian grids.
3 Numerical results
Figure 3: At left : geometry and boundary conditions of the ow over a circular cylinder. At right :
meshes from our simulations. One mesh line out of three is shown in both direction.
We have validated our method on both steady and unsteady ows around a circular cylinder in a
free-stream (see Fig. 3 (left)). The Reynolds number is based on the constant inlet velocity and
the cylinder diameter D. For all computations, the inlet boundary is located at X
u
= 8 length
units upstream of the obstacle ; the value of the time step is t = 0.01. We use non-uniform
Cartesian meshes rened in the wake of the cylinder as shown in Fig. 3 (right), with N
x
N
y
cells. In the vicinity of the cylinder, the cell distribution is uniform of size h = 0.04. For a
shape as simple as a cylinder prole with unit diameter, the level set function takes the simple
analytic expression :
(x, y) =
1
2

_
(x x
C
)
2
+ (y y
C
)
2
,
5
18
me
Congrs Franais de Mcanique Grenoble, 27-31 aot 2007
N
x
N
y
A X
d
L
w
u
min
C
D
S1 300 260 12 15 0.952 -0.032 2.025
S2 300 300 16 15 0.956 -0.032 2.017
Experiments [11] - -0.03 -
Linnick & Fasel [6] 0.93 - 2.06
Table 1: Summary of the steady computations at Re = 20 and comparison with results published in the
literature. L
w
: wake length, u
min
: velocity extremum in the near wake. C
D
: drag coefcient.
N
x
N
y
A X
d
St C
D
C
L
max
S1 300 260 12 15 0.169 1.341 0.350
S2 300 300 16 15 0.169 1.339 0.349
S3 336 300 16 20 0.169 1.339 0.349
Experiments [11] 0.16 0.17 1.21 1.41 -
Persillon & Braza [9] 0.165 1.253 0.395
Linnick & Fasel [6] 0.166 1.34 0.333
Table 2: Summary of the unsteady computations at Re = 100 and comparison with results published
in the literature. St : Strouhal number, C
D
: time average of drag coefcient. C
L
max
: peak value of
lift coefcient. All quantities have been computed from a time signal equal to 280 units, which gives a
frequency resolution equal to 1.8 10
3
.
where (x
C
, y
C
) is the cylinder center. This expression is discretized at cell corners prior to the
time integration.
For both steady ow at Re = 20 and time-periodic ow at Re = 100, we have studied the
inuence of the locations of the outow boundary X
d
and solid blockage A (see Tables 1 and
2). For Re = 100 in particular, it can be observed that grid independent results are reached for
A = 16 and X
d
= 15. These tables also compares our results with selected experimental and
numerical studies. Our computations show an overall good agreement with these results.
References
[1] O. Botella and Y. Cheny. On the treatment of complex geometries in a Cartesian grid ow solver with
the level set method. in European Conference on Computational Fluid Dynamics ECCOMAS CFD 06, P.
Wesseling, E. Oate & J. Priaux eds. (CD-ROM & WEB) 2006.
[2] O. Botella, M. Y. Forestier, R. Pasquetti, R. Peyret, and C. Sabbah. Chebyshev methods for the Navier-
Stokes equations : Algorithms and applications. Nonlinear Analysis, 47, 4157-4168, 2001.
[3] M. Drge and R. Verstappen. A new symmetry-preserving Cartesian-grid method for computing ow past
arbitrarily shaped objects. Int. J. Numer. Meth. Fluids, 47, 979-985, 2005.
[4] F. H. Harlow and J. E. Welch. Numerical calculation of time-dependent viscous incompressible ow of uid
with free surfaces. Phys. Fluids, 8, 2181-2189, 1965.
[5] S. Kang, G. Iaccarino, and P. Moin. Accurate and efcient immersed-boundary interpolations for viscous
ows. In Center for Turbulence Research Briefs, pp. 31-43, 2004.
[6] M. N. Linnick and H. F. Fasel. A high order immersed interface method for simulating unsteady incom-
pressible ows on irregular domains. J. Comput. Phys., 204, 157-192, 2005.
[7] R. Mittal and G. Iaccarino. Immersed boundary methods. Annu. Rev. Fluid Mech., 37, 239-261, 2005.
[8] S. Osher and R. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces. Springer, New-York, 2003.
[9] H. Persillon and M. Braza. Physical analysis of the transition to turbulence in the wake of a circular cylinder
by three dimensional Navier-Stokes simulation. J. Fluid Mech., 365, 23-88, 1998.
[10] R. W. C. P. Verstappen and A. E. P. Veldman. Symmetry-preserving discretization of turbulent ow. J.
Comput. Phys., 187, 343-368, 2003.
[11] M. M. Zdravkovich. Flow Around Circular Cylinders. Volume 1 : Fundamentals. Oxford University Press,
Oxford, 2003.
6