Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.

Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 762
The Numerical Simulation of Ship Waves Using Cartesian Grid
Methods
M.Sussman (Florida State University, USA)
D.Dommermuth (Science Applications International Corporation, USA)
ABSTRACT
Two different cartesian-grid methods are used to simulate the flow around the DDG 5415. The first technique uses a
“coupled level-set and volume-of-fluid” (CLS) technique to model the free-surface interface. The no-flux boundary
condition on the hull is imposed using a finite-volume technique. The second technique uses a level-set technique (LS) to
model the free-surface interface. A body-force technique is used to impose the hull boundary condition. The predictions
of both numerical techniques are compared to whisker-probe measurements of the DDG 5415. The level-set technique is
also used to investigate the breakup of a two-dimensional spray sheet.
1 INTRODUCTION
At moderate to high speed, the turbulent flow along the hull of a ship and behind the stern is characterized by
complex physical processes which involve breaking waves, air entrainment, free-surface turbulence, and the formation of
spray. Traditional numerical approaches to these problems, which use boundary-fitted grids, are difficult and time-
consuming to implement. Also, as waves steepen, boundary-fitted grids will break down unless ad hoc treatments are
implemented to prevent the waves from getting too steep. At the very least, a bridge is required between potential-flow
methods, which model limited physics, and more complex boundary-fitted grid methods, which incorporate more physics,
albeit with great effort and with limitations on the wave steepness. Cartesian-grid methods are a natural choice because
they allow more complex physics than potential-flow methods, and, unlike boundary-fitted methods, cartesian-grid
methods require minimal effort with no limitation on the wave steepness. Although cartesian-grid methods (CGM) are
presently incapable of resolving the hull boundary-layer, CGM can model wave breaking, free-surface turbulence, air
entrainment, spray-sheet formation, and complex interactions between the ship hull and the free surface, such as transom-
stern flows and tumblehome bows. The cartesian-grid methods that are described in this paper use the panelized geometry
that is used by potential-flow methods to automatically construct a representation of the hull. The hull representation is
then immersed inside a cartesian grid that used to track the interface. No additional gridding beyond what is already used
by potential-flow methods is required. We note that another variation of this approach is to use cartesian-grid methods to
track the free-surface interface and body-fitted grids to model the ship hull.
For the calculation of ship waves, VOF and level-set methods have certain advantages and disadvantages. VOF uses
the volume fraction (F) to track the interface. F=0 corresponds to gas and F=1 corresponds to liquid. For intermediate
values, between zero and one, there exists an interface between the gas and the liquid. The interface between the gas and
the liquid is sharp for a pure VOF method. Level-set methods use a level-set function to model the gas-liquid interface.
By definition, denotes gas, denotes liquid, and is the interface. For conventional level-set schemes, the
interface between the gas and the liquid is given a finite thickness[17], which is unlike conventional VOF schemes[3].
In the case of free-surface flows, where the density ratio between air and water is almost three orders
the authoritative version for attribution.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 763
of magnitude, the finite thickness of the interface that characterizes level-set methods has two advantages over VOF.
First, the finite thickness tends to smooth jumps in the tangential component of the velocity on the interface. Second, the
finite thickness tends to facilitate using multigrid methods to solve various types of elliptic equations that involve the
density.
The advection algorithm that is used for VOF conserves mass if the flow field is solenoidal. The level-set advection
equation tends to accumulate numerical errors. For the level-set method, the level-set function must be periodically
reinitialized to maintain a proper thickness for the interface, otherwise the interface would become either too thick or too
thin. The reinitialization process is a significant source of errors in the level-set method. Based on accuracy
considerations, the calculation of gravity-driven flows tends to favor VOF over level-set methods.
The interface is reconstructed from the volume fractions in VOF. During the reconstruction process, the interface
normals and curvature are calculated. Typically, the calculation of the interface normal and curvature are less accurate for
VOF than for level-set methods. The interface normals and curvature are calculated directly in level-set methods in terms
of gradients of the level-set function. As a result, the calculation of the normals and curvature are less costly for level-set
methods relative to VOF. The calculation of surface tension effects, which are a function of the curvature of the interface,
tends to favor level-set methods over VOF due to considerations of accuracy and efficiency.
On highly-stretched, multidimensional grids, VOF methods are less prone to aliasing errors than level-set methods.
Level-set methods incur errors as the interface rotates through highly resolved regions into regions that are not resolved
well. This type of aliasing error occurs in cartesian-grid methods when the mesh along one coordinate axis is more finely
resolved than along another coordinate axis.
By definition, the level-set and volume-of-fluid function both allow mixing of gas and liquid. This feature of level-
set and volume-of-fluid methods maybe desirable for modeling gas entrainment such as the air that is entrained by a
breaking wave. During the reinitialization process, level-set methods and “coupled level set and volume-of-fluid
methods” (CLS) use a signed distance function to update the level-set function and the thickness of the interface.
Naturally, the distance function could be used to model the intensity of turbulence and amount of gas entrainment as a
function of the distance to the interface.
Dommermuth, et al., (1998) used a stratified flow formulation to simulate breaking bow waves on the DDG 5415 at
a Froude number Fr=0.41. Their numerical results compared well to whisker-probe measurements in the bow region [8].
However, Dommermuth, et al., (1998) identified two issues that required further study. First, their stratified flow
formulation allowed the free-surface interface to become too diffuse. Second, the contact-line treatment didnot allow the
free surface to rise and fall cleanly along the side of the hull. The two new numerical approaches that are discussed in this
paper are attempts to remedy these problems.
Both numerical approaches use a signed distance function to represent the hull. The distance of a point to the hull is
negative inside the hull and positive outside the hull. The finite-volume approach uses the signed distance to calculate the
area and volume fractions for computational cells cut by the hull, whereas the body-force technique uses the signed
distance to prescribe a smooth forcing term. The coupled interface-tracking algorithm (CLS) uses level-set to calculate the
normals (and curvature if needed) to the free-surface interface that are used in VOF. The advection portion of the
algorithm is performed by VOF [16]. The level-set interface-tracking algorithm uses a new isosurface scheme to calculate
the zero level-set. Then the minimal distance between the cartesian points and the zero level-set is calculated in a narrow
band. The minimal distance is made positive in the water and negative in the air. This signed distance to the free surface
is used to reinitialize the thickness of the interface.
The two numerical approaches are used to simulate the flow around the DDG 5415. The CLS technique is still under
development, so only preliminary results are presented. The level-set technique includes upgrades to the numerical
technique that is described in [8]. Those upgrades include a new body-force formulation that is mollified, a new
reinitialization procedure, and a new finite-volume treatment of the convective terms. The original numerical procedure is
not mollified and does not use reinitialization. In addition, the original central-difference formulation of the convective
terms is not as robust as the new treatment using a flux integral formulation. We first review the governing equations and
then we discuss the numerical approaches. Finally, we present some preliminary numerical results which illustrate various
features of the numerical algorithms. The application of level-set methods to the breakup of spray sheets is also illustrated.
2 FIELD EQUATIONS
As in Dommermuth, et al., (1998), consider turbulent flow at the interface between air and water [8]. Let ui denote
the three-dimensional velocity field as a function of space (xi) and time (t). For an incompressible flow,
the authoritative version for attribution.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 764
the conservation of mass gives
(1)
ui and xi are normalized by Uo and Lo, which are the characteristic velocity and length scales of the body,
respectively. On the surface of the moving body (Sb), the fluid particles move with the body:
(2)
where Ui is the velocity of the body.
Let Vℓ and Vg respectively denote the liquid (water) and gas (air) volumes. Following a procedure that is similar to
for x ∈ Vg and for x ∈ Vℓ. The fluid interface
[13, 15], we let denote a level-set function. By definition,
corresponds to
The convection of is expressed as follows:
(3)
where d/dt=∂/∂t+ui∂/∂xi is a substantial derivative. Q is a sub-grid-scale flux which can model the entrainment of gas
into the liquid. Details are provided in [8].
Let ρℓ and µ ℓ respectively denote the density and dynamic viscosity of water. Similarly, ρg and µ g are the
corresponding properties of air. The flow in the water and air is governed by the Navier-Stokes equations:
(4)
where Re=ρℓUoLo/µℓ is the Reynolds number, is the Froude number, and is the Weber
number. g is the acceleration of gravity, and σ is the surface tension. Fi is a body force that is used to impose boundary
conditions on the surface of the body. P is the pressure. Ti accounts for surface-tension effects. δij is the Kronecker delta
symbol. As described in [8], τij is the subgrid-scale stress tensor. Sij is the deformation tensor:
(5)
ρ and µ are respectively the dimensionless variable densities and viscosities:
(6)
where λ=ρg/ρℓ and η=µg/µℓ are the density and viscosity ratios between air and water. For a sharp interface, with no
mixing of air and water, H is a step function. In practice, a mollified step function is used to provide a smooth transition
between air and water.
Based on [3, 4], the effects of surface tension are expressed as a singular source term in the Navier-Stokes equations:
(7)
where is the curvature of the air-water interface expressed in terms of the level-set function:
(8)
The pressure is reformulated to absorb the hydrostatic term:
(9)
where Pd is the dynamic pressure and Ph is a hydrostatic pressure term:
(10)
As discussed in [8], the divergence of the momentum equations (4) in combination with the conservation of mass (1)
provides a Poisson equation for the dynamic pressure:
(11)
where Σ is a source term. Equation 11 is used to project the velocity onto a solenoidal field.
3 ENFORCEMENT OF BODY BOUNDARY CONDITIONS
Two different cartesian-grid methods are used to simulate the flow around the DDG 5415. The first technique
imposes the no-flux boundary condition on the body using a finite-volume technique. The second technique
the authoritative version for attribution.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 765
imposes the no-flux boundary condition via an external force field. Both techniques use a signed distance function ψ to
represent the body. ψ is positive outside the body and negative inside the body. The magnitude of ψ is the minimal
distance between the position of ψ and the surface of the body.
With respect to the volume of fluid that is enclosed by the body (Vb), we define a function :
(12)
The function can be expressed in terms of a surface distribution of normal dipoles [10].
(13)
where n is the outward-pointing unit normal to the body, and R is a Rankine source, R=|x−x′|. ψ is expressed in terms
of as follows,
(14)
where |x−x′|min is the minimal distance between the field point x and the points on the body x′. In practice, the body
is discretized using triangular panels. As a result, the calculation of the minimal distance sweeps over all the triangles
comprising the body and must account for the possibility that the minimal distance may occur either at the corners of
triangle, along the edges of triangle, or inside the triangle.
3.1 Free-slip conditions
In the finite volume approach, the irregular boundary (i.e. ship hull) is represented in terms of ψ along with the
corresponding area fractions A and volume fractions V. V=1 for computational elements fully outside the body and V=0
for computational elements fully inside the body. The representation of irregular boundaries via area fractions and volume
fractions has been used previously in the following work for incompressible flows [1, 19, 5].
Recall the pressure equation,
(15)
with the following no-flow boundary condition:
(16)
where nwall is the outward normal drawn from the active flow region into the geometry region.
For each discrete computational element Ωi,j,k we define the geometry volume fraction V and area fraction A as
Γi+1/2,j,k represents the left face of a computational element; similar definitions apply to Γi−1/2,j,k, Γi,j+1/2,k,….
In order to discretely enforce the boundary conditions (16) at the geometry surface, we use a finite volume approach
for discretizing (15).
Given an irregular computational element Ωijk (see Figure 1), we have
The divergence theorem motivates the following second order approximation of the divergence ∇ · U at the centroid
of Ωijk:
(17)
In terms of geometry volume fractions Vijk and area fractions Ai+1/2,j,k, (17) becomes,
(18)
the authoritative version for attribution.
For a zero flux boundary condition at the wall, the last term in (18), is zero.
The finite volume approach, when applied to the divergence operator in (15) becomes:

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 766
and
Due to the no flow condition (16), the terms and cancel each other. The
resulting discretization for p is:
where, for example, (px)i+1/2,j,k is discretized as
3.2 No-slip conditions
The boundary condition on the body can also be imposed using an external force field. Based on Dommermuth, et
al., (1998), the distance function representation of the body (ψ) is used to construct a body force as follows:
Figure 1: Diagram of computational element (i, j) that is cut by the embedded geometry.
(19)
where cf is a friction coefficient. ∆ is used to mollify the body force such that it is gradually applied across the
surface of the body. Recall that ψ(x)≤0 corresponds to points within the body. Fi=0 outside of the body. is an
adjustment function:
(20)
To is the adjustment time. The adjustment function smoothly increases to unity from its initial value of zero. The
effect of the adjustment function is described in [8]. The adjustment function reduces the generation of nonphysical high-
frequency waves.
As constructed, the velocities of the points within the body are forced to zero. For a body that is fixed in a free
stream, this corresponds to imposing no-slip boundary conditions.
4 INTERFACE TRACKING
Two methods are presented in our work for computing ship flows. Both methods use a “front-capturing” type
procedure for representing the free surface separating the air and water. The first technique is based on the Coupled
the authoritative version for attribution.
volume-of-fluid and level set method (CLS) and the second technique is based on the level set method (LS) alone.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 767
4.1 CLS method
In this section, we describe the 2d coupled Level Set and Volume of Fluid (CLS) algorithm for representing the free
surface. For more details, e.g. axisymmetric and 3d implementations, see [16]. In the CLS algorithm, the position of the
interface is updated through the level set equation (level set function denoted by ) and volume of fluid equation (volume
fraction of liquid within each cell is denoted by Fij),
In order to implement the CLS algorithm, we are given a discretely divergence free velocity field uMAC defined on
the cell faces (MAC grid),
(21)
UMAC,
Given and we use a “coupled” second order conservative operator split advection scheme in order to
find and The 2d operator split algorithm for a general scalar s follows as
(22)
(23)
where Gi+1/2,j=si+1/2,jui+1/2,j denotes the flux of s across the right edge of the (i, j) th cell and
denotes the flux across the top edge of the (i, j)th cell. The operations (22) and (23) represent the case when one has the
“x-sweep” followed by the “y-sweep”. After every time step the order is reversed; “y-sweep” (done implicitly) followed
by the “x-sweep” (done explicitly).
The scalar flux si+1/2,j is computed differently depending on whether s represents the level set function or the
volume fraction F.
For the case when s represents the level set function we have the following representation for si+1/2,j (ui+1/2,j>0),
where
The above discretization is motivated by the second order predictor corrector method described in [2] and the
references therein.
For the case when s represents the volume fraction F we have the following representation for si+1/2,j (ui+1/2,j>0),
(24)
where
The integral in (24) is evaluated by finding the volume cut out of the region of integration by the line represented by
the zero level set of
The term found in (24) represents the linear reconstruction of the interface in cell (i, j). In other words,
has the form
(25)
A simple choice for the coefficients aij and bij is as follows,
(26)
(27)
The intercept cij is determined so that the line represented by the zero level set of (25) cuts out the same volume in
cell (i, j) as specified by In other words, the following equation is solved for cij,
the authoritative version for attribution.
where
and Fn+1 have been updated according to (22) and (23) we “couple” the level set function to the volume
After
fractions as a part of the level set reinitialization step. The level set reinitialization step replaces the current value of
with the exact distance to the VOF reconstructed interface. At the same time, the VOF reconstructed interface uses the
current value of to determine the slopes of the piecewise linear reconstructed interface.
Remarks:

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 768
• The distance is only needed in a tube of K cells wide K=ε/∆x+2, therefore, we can use “brute force” techniques
for finding the exact distance. See [16] for details.
• During the reinitialization step we truncate the volume fractions to be 0 or 1 if Although we truncate the
volume fractions, we still observe that mass is conserved to within a fraction of a percent for our test problems.
4.1.1 CLS Contact angle boundary conditions in general geometries
The CLS contact angle boundary conditions are enforced by extending into regions where Vij<1 (i.e. initializing
“ghost” values of in the inactive portion of the computational domain).
The contact angle boundary condition at solid walls is given by
(28)
where θ is a user defined contact angle and nwall is the outward normal drawn from the active flow region into the
geometry region.
In terms of (the free surface level set function) and ψ (the geometry level set function), (28) becomes
In figure 2, we show a diagram of how the contact angle θ is defined in terms of how the free surface intersects the
geometry surface.
Figure 2: Diagram of gas/liquid interface meeting at the solid. The dashed line represents the imaginary interface created thru
the level-set extension procedure.
The “extension” equation has the form of an advection equation:
(29)
In regions where is left unchanged.
For a 90 degree contact angle (the default for our computations), we have
In other words, information propagates normal to the geometry surface.
For contact angles different from 90 degrees, the following procedure is taken to find uextend:
Remarks:
• In 3d, the contact line (CL) is the 2d curve which represents the intersection of the free surface with the geometry
surface (ship hull). The vector n2 is orthogonal to the contact line (CL) and lies in the tangent plane of the
geometry surface.
and ψ are defined within a narrow band of the zero level set of we can also define uextend within
• Since both
a narrow band of the free surface.
• We use a first order upwind procedure for solving (29). The direction of upwinding is determined from the
extension velocity uextend. We solve (29) for τ=0…ε.
the authoritative version for attribution.
• For viscous flows, there is a conflict between the no-slip condition and the idea of a moving contact line. See [6,
8, 9, 12] and the references therein for a discussion of this issue. We have performed numerical studies for
axisymmetric oil spreading in water under ice [18] with good agreement with experiments. In the future, we wish
to experiment with appropriate slip-boundary conditions near the contact line.
4.2 Level-set method
A key part of level-set methods is reinitialization. Without reinitialization, the thickness of the interface between the
gas and the liquid can get either too thick or

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 769
too thin. Reinitialization is based on the construction of a signed distance function that represents the distance of points
from the gas-liquid interface. By definition, the signed distance is positive in the liquid and negative in the gas. At the
interface, the distance function is zero. A variety of methods have been utilized for calculating the signed distance
function, including a hyperbolic equation [17] and direct methods [16]. The hyperbolic equation methods tend to be less
accurate but more efficient than direct methods. Here, we outline a direct method that can be efficiently implemented on
parallel computers with second-order accuracy. The numerical scheme can also be generalized to higher order.
First, calculate the intersection points (xp) where the zero level-set crosses each of the cartesian axes. At these
intersection points calculate the normal to the interface (np). Together, xp and np determine local approximations to the
planes that pass through the zero level-set. For points that are within a narrow band of these planes, calculate the minimal
distance to the planes. Once the minimal distance is calculated, assign the sign of the distance function based on the sign
of the level-set function.
For example, consider a zero crossing along the z-axis. Locally, near the zero crossing, the level-set function is fitted
with Lagrange polynomials.
(30)
where zo is offset where interpolated level-set function Lk are Lagrange polynomials and are discrete values
of near the zero level-set along the z-axis. K−1 is the degree of the interpolating polynomial. zo is calculated directly
for low-order polynomials and iteratively for high-order polynomials. Let xo=(xo, yo, zo), where (xo, yo, zo) is the coordinate
of the zero crossing.
The unit normal no at the zero crossing is calculated in terms of the level-set function:
(31)
where the gradient terms are calculated using finite difference formulas of desired order.
The minimal distance (s) between a point (xp) and a plane lies along the unit normal to the plane. Denote the position
where the point intersection occurs as xs, then
(32)
where s is expressed in terms of a dot product:
(33)
Note that higher-order corrections involve curvature terms, etc. As long as |xo−xs|≤∆g, where ∆g is the grid size, then
s is potentially the minimal distance to the zero level-set. Other candidates include planes in the neighborhood of xp. On a
structured grid, shifts along the cartesian axes can be performed to consider other candidates. Only xp near the zero level-
set are required in the reinitialization procedure. A simple procedure for finding points near the zero level-set involves
weighted averages. First construct a stair-case approximation (Φ) to the zero level-set:
(34)
A weighted average along the k-th indice is
(35)
Similar expressions hold along the i-th and j-th indices. Repeated applications of weighted averages provide a
narrow band that encompasses the zero level-set. The narrow band corresponds to the region The signed
distance function D is expressed in terms of the level-set function and the minimal distance:
(36)
Based on [17], is reinitialized as follows:
(37)
where ∆ is the desired thickness of the interface.
5 FLUX INTEGRAL METHODS
We define the temporal and spatial averaging over a time step and a cell as follows:
(38)
where here, the tilde and overbar symbols respectively denote temporal and spatial averaging. ∆t is the time step, and
∆V is the volume of the cell.
the authoritative version for attribution.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 770
As an example, consider the application of the preceding operator to the level-set equation (3):
(39)
where here superscript n denotes the time level. We focus our attention on the convective term. The convective term
accounts for the flux of the level-set function across the faces of the control volume. A second-order approximation for
the flux across one face of a cell is provided below:
(40)
where is the flux across the positive face along the x axis. The limits of integration are provided below:
(41)
where ∆x, ∆y, and ∆z are the lengths of the cell along the cartesian axes. u+ is the normal component of fluid velocity
at the center of positive face along the x-axis. v+ and v− are the normal components of the fluid velocities at the centers of
the positive and negative faces along the y-axis. Similar definitions hold for w+ and w−.
In a mapped coordinate system, the expression for the flux is
(42)
where J is the jacobian, and x, y, and z are functions of r, s, and t:
(43)
For this particular approximation, the jacobian is
(44)
On any one face the stencil associated with the Lagrangian interpolation of is 3×3×3=27 points, but for the
entire cell, the stencil is 5×5×5=125 points. We use a upwind-biased stencil for the momentum equations and a symmetric
stencil for the level-set function. The diagonal and cross terms in the momentum equations are treated the same.
Generally, we use eight-point Gaussian quadrature to evaluate the flux over each face. Details of the numerical algorithm
are described in [7]. Various types of limiters are described in [11].
6 PRELIMINARY RESULTS
In section 6.1, we present preliminary computations of flow past a DDG 5415 ship. In section 6.2, we present
preliminary computations of the breakup of a two-dimensional spray sheet.
6.1 Ship Wave Results
As a demonstration of the level-set and the coupled level-set and volume-of-fluid formulations, we predict the free-
surface disturbance near the bow of the DDG 5415 moving with forward speed. The experiments were performed at the
David Taylor Model Basin (DTMB), and are available via the world wide web at http://www50.dt.navy.mil/5415/. This is
the same flow that Dommermuth, et al., (1998) originally investigated using their stratified flow formulation [8]. As
before, we only consider the high speed case. For this case, a plunging breaker forms near the bow. Air is entrained and
the authoritative version for attribution.
splash up occurs where the wave reenters the free surface. There is flow separation at the stern, and the transom is dry. A
large rooster tail forms just behind the stern.
Based on the speed (Uo=6.02Knots) and the length (Lo=5.72m) of the model, the Reynolds and Froude numbers are
Re=1.8×107 and Fr2=0.41. The effects of surface tension are not included. The density ratio of air and water is λ=0.0012
and the ratio of the dynamic viscosities is η=0.018.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 771
In regard to the numerical parameters for the level-set formulation, we use a friction coefficient cf=500 in the body-
force term (19). The adjustment time is To=0.02. For the level-set formulation, the length and width of the computational
domain are L=2.5 and W=1.50. The height of the air above the mean free-surface is h=0.15 and the depth below the mean
free-surface is d=1.0. One grid resolution is used with 512×128×129 grid points. Three different levels of grid stretching
are used along the y- and z-axes. For the highest grid resolution, the smallest grid spacing is 2.6×10−3 along the y-axis and
3.6×10−4 along the z-axis. For the medium resolution simulation, the smallest grid spacing is 3.8×10−3 along the y-axis
and 1.8×10−3 along the z-axis. For the coarsest grid simulation, the smallest grid spacing is 3.8×10−3 along the y-axis and
3.5×10−3 along the z-axis. The grid spacing (4.9×10−3) is constant along the x-axis for all three cases. The thicknesses of
the free-surface interfaces for the fine, medium, and coarse simulation are respectively ∆=0.05, 0.025, and 0.0125. The
durations of the coarse and medium resolution simulations are t=0.76 and t=0.68, respectively. No special treatment is
used for the level-set function inside the ship. These durations correspond to about three quarters of a ship length based on
the present normalization. For these durations, the flow is steady near the bow and still evolving near the stern. (The fine
resolution simulation is still evolving, and it is not possible at this time to present complete results. More complete results
will be provided at the symposium and in the discussion section of this paper1.) The ship is centered in the computational
domain with the same fixed sinkage and trim as used in the experiments. In order to construct the body force term, the
hull is panelized using approximately 4000 panels.
Coarse and medium resolution simulations have been performed using the CLS formulation. The coarse simulation
uses 256×64×64 grid points, and the fine resolution uses 512×128×128 grid points. The length, width, and height of the
computational domain are L=2, W=0.5, and H=0.5, respectively. The water depth is d=0.25. The grid spacing is constant
along all three cartesian axes. In the next phase of our research, we will implement grid stretching, which will allow
greater water depths to be simulated. The durations of the CLS simulations are t=0.75. Unlike the level-set results, the
CLS results extend the free-surface interface into the hull using the techniques outlined earlier in our paper.
The free-surface elevation was measured at DTMB using a whisker probe. Twenty-one transverse cuts were
performed near the bow, extending from x=0 to x=0.178 in dimensionless units. The whisker probe measures the highest
point of the free surface. In regions where there is wave breaking, the whisker probe measures the top of the breaking
wave. Seventeen transverse cuts were performed in the stern, extending from x=1.01 to x=1.22.
Figures 3 and 4 compare measurements at the bow and stern to the numerical predictions. The bow measurements
include profile and whisker-probe measurements. Comparisons to the bow data are performed at four stations: x=0.0444,
x=0.0622, x=0.0800, and x=0.0978. The circular symbol denotes profile measurements. The solid black lines denote the
outline of the hull and the whisker-probe measurements. The solid blue line is medium CLS and the dashed blue line is
coarse CLS. The solid red line is medium level-set and the dashed red line is coarse level-set. In general, the CLS
technique captures the rapid rise up the side of the hull. The level-set technique does less well in this regard. In the outer-
flow region the CLS coarse results are slightly better than the CLS fine results. This may be attributed to the shallow
depth that is used in the CLS. The level-set results appear to converge better in the outer-flow region, but the results of the
fine simulation are required for confirmation.
Figure 4 shows the entire flow around the ship for the medium resolution level-set simulation. The stern whisker-
probe measurements are overlaid for the purposes of comparison. Although the numerical results are not stationary, the
shape of the stern contours show general agreement with laboratory measurements. However, the amplitude of the
numerical results are significantly lower than the measurements. Note that the stern is partially dry in the numerical
simulations. The outline of the hull is visible in the numerical simulations because the level-set function intersects the hull.
6.2 Spray Sheet Results
The Navier-Stokes equations in combination with a level-set formulation are used to study the breakup of two-
dimensional sheet of water. The sheet is lo=6mm thick. The length of the sheet is 24mm. The top and bottom of the sheet
are bounded by air. The initial mean-velocity of the water is uo=3m/s. The initial rms turbulent velocity of the water is
The air is initially quiescent. Based on the sheet thickness (lo) and the mean velocity (uo), the Reynolds
number is Re=uolo/µ=18,000 and the Weber number is where µ is the kinematic viscosity of water, ρ
is the water density, and σ is the surface tension. The density and viscosity ratios
the authoritative version for attribution.
1We would have performed longer simulations, but the NAVO T3E was unexpectedly shutdown for five days of maintenance just
before this paper was due.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 772
are λ=ρg/ρℓ=0.0012 and η=µg/µℓ=0.018, which are appropriate for air-water interfaces. This parameter regime roughly
corresponds to experiments that were performed by Sarpkaya and Merrill (1998), [14]. Numerical convergence is
established using 20482 and 40962 grid points. Second-order accuracy in space is established. A third-order Runge-Kutta
scheme is used to integrate the system of equations with respect to time. Mass is conserved to within 0.25% throughout
the entire calculation.
Figure 5 illustrates the evolution of a two-dimensional spray sheet. The black contour lines indicate the interface
between air and water. The water sheet is bounded by air both at the top and the bottom of the sheet. The color contours
denote the vorticity. The flow is turbulent within the water sheet and laminar in the air. The mean velocity and rms
velocity profiles are initially top-hat functions. The flow is moving from left to right. The turbulent fluctuations in the
water are initially immersed below the top of the sheet and above the bottom of the sheet (see Fig. 5: t=0).
The turbulence in the water diffuses and interacts with the interfaces (see Fig. 5: t=2.5). The initial interaction is a
roughening of the air-water interface. A thin boundary layer forms in the air. The boundary layer is colored blue
(negative) at the top of the sheet and colored red (positive) at the bottom of the sheet. As the interface gets rougher and
ligaments begin to form, the air separates from the back of the ligaments. The boundary layer thickens, and air is dragged
along the top and the bottoms of the sheet.
Primary vortex shedding initially occurs behind the ligaments (see lower left of sheet in Fig. 5: t=5). As the primary
vortices are shed, their interactions lead to the formation of secondary and tertiary vorticity (see upper middle of sheet in
Fig. 5: t=7.5). Vortices are periodically shed from the backs of ligaments (see lower middle of sheet in Fig. 5: t=10).
There is evidence of vortex merging both in the air and in the water (see upper left of Fig. 5: t=17.5). Although there is
significant flow separation in the air, there is little or no separation in the water. The largest ligaments are formed by
eddies impinging on the interface (see upper left of Fig. 5: t=12.5). Cavities form in regions where primary vortices are
trapped. The inlets to the cavities shed secondary vorticity, which tends to make the cavities even larger (see middle of
sheet in Fig. 5: t=15). At the inlets to the cavities, vortex pairs are formed. Under their own self-induced velocities, the
vortex pairs move into the cavities where they diffuse.
Note that droplets do not actually form at the tips of the ligaments because 2d flows are not subject to the same
instabilities as 3d flows. The turbulent kinetic energy tends to concentrate in the thicker portions of the deformed spray
sheet. The flow within the ligaments is relatively benign. In agreement with theory, the pressure at the tips of the longest
ligaments roughly scales like P=(Wer)−1 where r is the radius of curvature of the tip.
7 CONCLUSION
In this paper, we have outlined the key numerical algorithms for simulating free-surface flows on cartesian grids
using level-set and coupled level-set and volume-of-fluid techniques. Preliminary numerical results have been shown for
ship waves and spray sheets. The ship wave results indicate that cartesian-grid methods are capable of resolving the flow
around a ship if the grid resolution is sufficient. Near the bow and stern, we estimate that the grid spacing along all three
cartesian axes should be ∆=0.0005 (based on ship length) in order to resolve breaking waves. On a parallel computer, it is
possible to approach this level of grid resolution, but adaptive gridding may also be required to fully resolve the entire
flow around a ship [15]. Alternatively, cartesian-grid methods could be embedded in more conventional boundary-fitted
methods to capture complex flows near the bow or stern. The spray-sheet results show that cartesian-grid methods are
capable of resolving the air and water boundary layer at realistic Reynolds numbers.
Acknowledgments. The first author is supported in part by NSF Division of Mathematical Sciences under award
number DMS 9996349. The second author is supported by ONR under contract number N00014–97-C-0345. Dr. Edwin
P.Rood is the program manager. The numerical simulations have been performed on the T3E computer at the Naval
Oceanographic Office using funding provided by a Department of Defense Challenge Project. We are very grateful to Mr.
George Innis, Dr. James Rottman, and Mr. Andrew Talcott for assistance with this paper.
REFERENCES
[1] A.S.Almgren, J.B.Bell, P.Colella, and T.Marthaler. A cartesian grid projection method for the incompressible euler equations in complex geometries.
SIAM J. Sci. Comput., 18(5):1289–1309, 1997.
[2] J.B.Bell, P.Colella, and H.M.Glaz. A second-order projection method for the incompressible Navier-Stokes equations. J. Comput. Phys., 85:257–
283, December 1989.
the authoritative version for attribution.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 773
[3] J.U.Brackbill, D.B.Kothe, and C.Zemach. A continuum method for modeling surface tension. J. Comput. Phys., 100:335–353, 1992.
[4] Y.C.Chang, T.Y.Hou, B.Merriman, and S.Osher. Eulerian capturing methods based on a level set formulation for incompressible fluid interfaces. J.
Comput. Phys., 124:449–464, 1996.
[5] P.Colella, D.T.Graves, D.Modiano, E.G.Puckett, and M.Sussman. An embedded boundary/volume of fluid method for free surface flows in irregular
geometries. In proceedings of the 3rd ASME/JSME joint fluids engineering conference, number FEDSM99–7108, San Francisco, CA, 1999.
[6] R.G.Cox. The dynamics of the spreading of liquids on a solid surface, part 1. viscous flow. J. Fluid Mech., 168:169–194, 1986.
[7] D.G.Dommermuth. Numerical Flow Analysis (nfa) working papers. Technical report, Science Applications International Corporation, 2000.
[8] D.G.Dommermuth, G.E.Innis, T.Luth, E.A. Novikov, E.Schlageter, and J.C.Talcott. Numerical simulation of bow waves. In proceedings of the
Twenty Second Symposium on Naval Hydro., pages 508–521, Washington, D.C., 1998.
[9] L.M.Hocking and A.D.Rivers. The spreading of a drop by capillary action. J. Fluid Mech., 121:425–442, 1982.
[10] H.Lamb. Hydrodynamics. Dover Publications, New York, 1932.
[11] B.J.Leonard. Bounded higher-order upwind multidimensional finite-volume convection-diffusion algorithms. In W.J.Minkowycz and E.M.Sparrow,
editors, Advances in Numerical Heat Transfer, volume I, pages 1–58. Taylor & Francis, 1997.
[12] C.G.Ngan and E.B.Dussan V. On the dynamics of liquid spreading on solid surfaces. J. Fluid Mech., 209:191–226, 1989.
[13] S.Osher and J.A.Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations. J. Comput. Phys.,
79(1):12–49, 1988.
[14] T.Sarpkaya and C.Merrill. Spray formation at the free surface of liquid wall jets. In proceedings of the Twenty Second Symposium on Naval Hydro.,
pages 796–808, Washington, D.C., 1998.
[15] M.Sussman, A.Almgren, J.Bell, P.Colella, L.Howell, and M.Welcome. An adaptive level set approach for incompressible two-phase flows. J.
Comput. Phys., 148:81–124, 1999.
[16] M.Sussman and E.G.Puckett. A coupled level set and volume of fluid method for computing 3d and axisymmetric incompressible two-phase flows.
J. Comp. Phys. accepted for publication.
[17] M.Sussman, P.Smereka, and S.J.Osher. A level set approach for computing solutions to incompressible two-phase flow. J. Comput. Phys., 114:146–
159, 1994.
[18] M.Sussman and S.Uto. Computing oil spreading underneath a sheet of ice. Technical Report CAM Report 98–32, University of California, Los
Angles, July 1998.
[19] H.S.Udaykumar, H.C Kan, W.Shyy, and R.Tran-Son-Tay. Multiphase dynamics in arbitrary geometries on fixed cartesian grids. J. Comput. Phys.,
137(2):366–405, 1997.
the authoritative version for attribution.

OCR for page 762

About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
the authoritative version for attribution.
Figure 3: Flow near bow.
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS
774

OCR for page 762

About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
the authoritative version for attribution.
Figure 4: Flow near stern.
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS
775

OCR for page 762

About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
the authoritative version for attribution.
Figure 5:2d spray sheet.
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS
776

OCR for page 762

About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
the authoritative version for attribution.
Figure 5:2d spray sheet continued.
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS
777

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 778
DISCUSSION
U.Bulgarelli
Instituto Nazionale per Studi ed Esperienze di
Architettura Navale, Italy
In your algorithm do you have already adopted the adaptive grid in the 3D geometry?
AUTHOR'S REPLY
We are in the process of developing a body-fitted method. The method will be described at the next symposium.
DISCUSSION
K.Hendrickson
Massachusetts Institute of Technology, USA
The authors of this paper show an aggressive use of a numerical method which has, to date only, been used for smaller
engineering problems. The fact that they are attempting it for this type of problem says much about their patience and
ambition. The blending of the volume of fluid and level set methods is quite creative and shows promising results. I believe
that the Cartesian Grid Method is wonderfully useful in that many of the gridding difficulties have been removed and/or
reduced to panel method that has been dealt with in detail in the literature. I wish them both luck as they push the method
further.
QUESTIONS
1. At this stage in development of the two methods, it seems that the coupled level-set/volume of fluid
technique (CLS) is more accurate/robust in treating the hull boundary conditions mainly because it can better
define where the hull lies in the Cartesian grid. Have the authors done any investigation on the effects of the
mollified body force term used in the level-set (LS) technique? In using comparisons to the waterline
measurements as the benchmark, the exact hull position would likely be a critical point. Is it possible that the
mollified body force term is smoothing out the hull to the extent that it is affecting the waterline results?
Would less mollification or higher resolution in the region near the body produce better LS results? This can
almost be inferred from Figure 3 in the paper.
2. The choice of friction coefficient in the body force term (equation 19) seems to be somewhat arbitrary. Have
the authors done any type of parametric study on a range of friction coefficients and their effect on the LS
results?
3. Considerable effort has been invested in the LS community to address reinitialization, which is also done in
this paper. The reinitialization issue comes about because a Lagrangian thought process has been applied to a
Eulerian method. In most LS formulations, the advection of the level-set function, ~, is performed using the
velocity of the fluid. This causes the LS function to lose its distance function property and require
reinitialization. It is possible to construct a velocity field such that the distance function remains one [1].
Have the authors considered this type of LS formulation?
4. What is the computational cost comparison between the CLS and LS methods at the resolutions submitted in
the paper?
5. How do the authors feel CGM compare to other less computationally expensive capabilities such as unRANS
or 2D+T methods?
6. What do the authors consider to be the major limitations of the CGM, both CLS and LS, in terms of their
applicability to Marine Hydrodynamics and Computational Ship Hydrodynamics?
REFERENCES
1. Adalsteinsson, D. and Sethian, J.A. “The Fast Construction of Extension Velocities in Level Set Methods,” J. Comp. Physics, Vol. 148, 1999, pp. 222.
AUTHOR'S REPLY
1. The coupled level-set/volume (CLS) of fluid technique has a more accurate treatment of the hull boundary
condition than the level-set method (LS). Our tests indicate that mollified body-force terms improve convergence.
2. The body-force term is as large as possible without violating the Courant condition.
3. The advantage of our reinitialization procedure is its accuracy, which can be generalized to any order. Other
procedures,
the authoritative version for attribution.

OCR for page 762

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
THE NUMERICAL SIMULATION OF SHIP WAVES USING CARTESIAN GRID METHODS 779
such as the one proposed by the discussor, are effective away from the interface.
4. The CLS method is about twice as expensive as the LS method. However, the computational costs associated with
both methods are less than ten percent of the Poisson solver.
5. Interface tracking methods are capable of capturing physics that unRANS and 2D+T will never be capable of
modeling. Although interface tracking methods are more computationally expensive than unRANS and 2D+T, this will
become less of an issue as computers become faster. Ten years from today, interface tracking will be the method of
choice for modeling breaking waves and the near-field flow around naval combatants.
6. The treatment of the hull boundary condition is not accurate enough. This issue is currently being addressed by
using a body-fitted grid with a level-set treatment of the free-surface elevation in the near field of the hull. In the outer-
flow region, the inner solution is matched to a combination of spectral methods and panel methods. This matching
procedure reduces the number of grid points and the amount time that is required to generate 3D grids.
the authoritative version for attribution.