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 833
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 833
Simulation of UUV Recovery Hydrodynamics
S.Huyer, J.Grant (Naval Undersea Warfare Center, USA)
ABSTRACT:
A novel method to compute the 3-D unsteady hydrodynamics with application to undersea vehicles is presented.
This approach solves the vorticity equation, which is derived from the momentum equation of the Navier-Stokes
equations. Most problems of Navy interest involve incompressible flow, which can be described in terms of the vorticity
alone. Velocity is an integral quantity of the instantaneous vorticity field. Specific geometries are represented using
surface source and vortex panels whose strength is prescribed to satisfy the no-slip and no-flux boundary conditions.
Vorticity is diffused from the vortex sheets onto the body surface to maintain a vorticity balance. Vorticity in the flow is
specified at points and the vorticity at any other point in the field is obtained via linear interpolation. Interpolation is
performed by constructing tetrahedra using Delaunay triangularization. Tetrahedra provide the control volume to integrate
over to obtain the velocity and the connectivity of the control points provides a basis to construct derivatives. A Baldwin-
Lomax eddy viscosity model was implemented into the solution algorithm to model turbulent flow effects. This method
was then validated for two disparate flow cases—flow past an unmanned undersea vehicle (UUV) at Reynolds numbers
of one million and unsteady flow development past a cone. Attached flow past the UUV was compared with empirical
turbulent flat plate results. Quality of the flow past a cone was compared with data obtained with experimental data.
Validation of this method allows for a subsequent simulation of a UUV recovery problem.
INTRODUCTION:
Undersea vehicle hydrodynamics pose significant challenges for the computation of complex, three-dimensional
unsteady flow fields. Major examples include submarine maneuvering problems, low speed maneuvering and control and
unmanned undersea vehicle (UUV) recovery. Many of these complex flows are characterized by the production of
vorticity and its subsequent interaction with the vehicle. This is generally the case when a vehicle encounters the wake of
another body. For cases where vorticity is dominant, a Lagrangian vorticity based method can be used to compute the
complex unsteady two-body hydrodynamics. Here, incompressible, unsteady fluid flow can be characterized by the
instantaneous vorticity field alone. For cases of flow past bodies, the vorticity distribution in the boundary layer and wake
determines the characteristics of the entire flow field. Vorticity based methods would seem to be a natural fit to solving
these types of problem. This technique utilizes the velocity-vorticity formulation of the Navier-Stokes equations to solve
for these flow variables on a Lagrangian mesh (calculational points advected by the local flow). The solution
methodology has distinct advantages over traditional methods that rely on fixed grid solutions. The Lagrangian vorticity
method is essentially grid free. It does not rely on a grid (structured or unstructured) in the traditional sense to evolve the
vorticity associated with the unsteady flow. Vorticity is continually generated at the surface and is represented by the
computational points that are also continuously generated. This makes the method naturally adaptive to coherent vortex
structures and since the vorticity is advected by the flow, it is subject to little numerical diffusion. The use of a
Lagrangian mesh allows for the straightforward treatment of moving surfaces and does not require incorporation of
additional terms due to non-inertial reference frames. In addition, multiple bodies are included in a straightforward
manner. The traditional vortex method (Chorin, 1973) describes the vorticity field by means of isotropic elements or
‘blobs', which have a strength that depends only on distance from their center; a frequently used strength distribution is
the Gaussian. Careful comparison with theoretical and experimental data of two-dimensional calculations using isotropic
blobs of uniform size has been reported by Sethian and Ghoniem (2) for a backward-facing step. High-resolution 2-D
computations for flow past a cylinder have been conducted by Koumoutsakos and Leonard (1995) and Subramaniam
(1996) and represent the standard for using blob methods to compute unsteady flow past surfaces. The original method
has the feature that the identity and location of neighboring blobs are not needed to compute the Biot-Savart integral
(which determines the velocity), so the algorithm for this computation is simple. Later works (e.g., Koumoutsakos and
Leonard, 1995) typically employ accelerated methods (Greengard and Rokhlin, 1987; Strickland and Baty, 1995) to avoid
an order N2 calculation (where N is the number of elements). Even with these elaborations, the simplicity of an approach
the authoritative version for attribution.
based on the traditional method remains attractive.
There remain significant difficulties when applying vortex blobs to flows past a surface, however. One is the fact
that near the surface, the blob vorticity distribution (a Gaussian) actually penetrates the surface so a finite value of
vorticity is on the inside of the surface. This can cause problems when computing the boundary conditions. Another is
that without sufficient overlap of the blob radii, the computed velocity and associated vorticity field will be extremely
noisy. In the current method, which will be presented later, the vorticity in the field is piecewise
OCR for page 834
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 834
continuous and is linearly interpolated between points in the field. No overlap of elements is necessary to maintain a
smooth velocity and vorticity field. In contrast, blob methods rely on a superposition of the blob functions to determine
the field vorticity and subsequent velocity distribution. If sufficient overlap of the functions is not maintained, a noisy
velocity distribution can result causing poor satisfaction of the boundary conditions and a diverging solution. This overlap
can be very difficult to maintain, especially near the surface as the blobs are advected by the flow. For this reason, most
blob methods (e.g. Koumoutsakos and Leonard, 1995) employ periodic interpolation of the blob strength onto a regular
grid to most efficiently maintain overlap to avoid noisy boundary conditions. This can lead to artificial numerical
diffusion, however, and still does not solve the problem of vorticity penetration into the surface due to the blob function.
Another difficulty with blob methods is the wide range of scales that emerge during the evolution of flows past a
surface. An example is the anisotropy of the vorticity distribution in a thin boundary layer. For these cases, the gradient of
vorticity in the normal direction is much greater than that in the tangential direction. This presents another compelling
motivation for developing an algorithm not based on elements of uniform size, but one in which the size and shape of the
elements adapt to the local spatial distribution of the vorticity field. Blobs, whatever their size, have a constant radius and
are therefore isotropic in nature. For thin boundary layers, anisotropic elements are desired. As problems are extended to
three-dimensions, the above problems are only exacerbated. For this reason, there are only minimal examples of vortex
blob solutions for flow past surfaces (e.g. Gharakhani and Ghoniem, 1996)
We have developed a novel technique directly solving the vorticity equation on a Lagrangian mesh. Vorticity is
specified at points in the field and is linearly interpolated by constructing tetrahedra. This method has the advantage that
the elements are connected allowing the vorticity to be locally approximated as a function of position. This approach was
introduced by Russo and Strain (1994) who examined inviscid vorticity fields for two-dimensional flow. Huyer and Grant
(2000) extended the method to examine viscous flow past streamlined and bluff bodies. The current approach treats fully
three-dimensional flow past multiple bodies. The vorticity is determined on a number of points in the field. These points
are then connected to form a set of tetrahedra via a Delaunay tetrahedralization algorithm. The vorticity at any other point
in the region is then approximated by linearly interpolating the vorticity at the nodal points with given shape functions
based on the geometry of the tetrahedra. First and second order derivatives are computed via a second order least squares
formulation from the elements connected at a given node (Marshall and Grant, 1995). The vorticity field and the locations
of the calculation points are updated at each time step. Since the formulation is a Lagrangian approach, the advection term
is automatically included. Viscous diffusion is accomplished using both an effective diffusion velocity (i.e. Strickland and
Baty, 1995) as well as the second order Laplacian. Integration of the Biot-Savart integral provides the velocities.
Another difficulty in treating high Reynolds number flow is accounting for turbulence. Traditional vortex methods
have relied on random walk and other stochastic methods to simulate the high Reynolds number turbulent flow. Thus far,
most deterministic solutions for viscous flow have been limited to low Reynolds number laminar cases. In the present
method, a deterministic approach utilizing the full viscous equations was desired. Therefore, a turbulence model has been
introduced to account for the lack of spatial resolution required to properly treat fully turbulent flow directly. A Baldwin-
Lomax turbulence model was implemented into the viscous solution methodology to better model the turbulent flow
characteristics. To the author's knowledge, this has yet to be implemented in any vorticity based solution methodology.
This paper includes unsteady computational data collected for flow past a UUV at a Reynolds number of one million
and unsteady flow past a cone for a Reynolds number of 50,000. Attached turbulent boundary layer flow over the UUV is
compared with turbulent velocity profiles obtained from empirical results by Spalding and Coles (see White, 1974). Time
averaged wake velocity data for flow past a cone is compared with data presented by Calvert (1967). The docking cone
and UUV are fixed in space and the flow field computed. These validation tests are performed to establish a high
confidence level in the method. This is followed by a simulation of UUV recovery and discussion of the results.
METHODOLOGY:
Surface Definition:
An example of an Unmanned Undersea Vehicle (UUV) surface mesh is shown in Figure 1. To construct the mesh,
surface body points and unit normals are required. Pseudo points are placed just above the body points. A Delaunay
tetrahedralization algorithm (Borouchaki and Lo, 1995) is then applied. This essentially constructs tetrahedra
encompassing all points and leaves a convex surface. The surface then consists of the faces of the tetrahedra, which are
exclusively made up of the original body points. The UUV surface was defined with 586 points and a total of 1150
surface panels. These panels are then used to define the surface source and vortex panels.
Satisfaction of the Surface Boundary Condition:
Each panel on the body surface carries two velocity generators: a surface vortex distribution lying in the plane of the
panel and a potential source. The sources are needed, mathematically, to ensure the no-flux boundary condition is met
the authoritative version for attribution.
properly. Uhlman and Grant (1993) showed that as the number of panels increases to infinity, the source strength
approaches zero since the surface vorticity distribution can satisfy the no-slip and no-flux boundary conditions
simultaneously. Both distributions are taken to be uniform over an individual panel and lie in an infinitely thin sheet on
the surface. Thus the vortex strength parameter characterizing a panel is the velocity jump across the panel. The velocity
due to a potential source, σ, and vortex panel strength, on a surface S is:
OCR for page 835
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 835
(1)
(2)
Figure 1: Unmanned Undersea Vehicle
Geometry
where σn and are the discretized strengths for a panel n of surface are Sn and:
The no-slip and no-flux boundary conditions are then applied to the system of equations to solve for the source and
vortex sheet strengths.
Initial Volume Vorticity Distribution:
After the surface panels are defined, it is required to initialize the volume vorticity beginning with the surface
vorticity. Nodal vorticity values are located at the surface body points used to define the geometry. Additional points are
then placed in layers staggered over the body nodes and the surface panel centroids. These layers are located at a constant
normal distance initially and are separated by a distance according to a number of boundary layer thicknesses. A single
thickness is defined based on a viscous diffusion length scale of where ν is the kinematic viscosity and ∆t the time
step. The initial layer is very thin and is of the order of 10% of a flat plate boundary layer thickness
Linearized Tetrahedral Vorticity Elements:
In the present method, nodal vorticity values are known and a linear variation of vorticity between nodal points is
assumed. At each time step, a Delaunay triangularization routine is used to form an unstructured mesh connecting each
nodal point thus forming tetrahedral elements. Delaunay triangulation effectively optimizes the aspect ratio of all
tetrahedra constructed from a random distribution of points. For more detail concerning this method, the reader is referred
to Borouchaki and Lo (1995). For a single tetrahedral element, Zienkiewicz (1977) derives the four shape functions that
are described as a function of the geometric location of the four vertices. The shape function values are 1.0 at their
respective nodes and 0.0 at each of the other three nodes. The vorticity over the element can then be expressed as:
(3)
The velocity (from the Biot-Savart Integral) is:
(4)
In order to solve this integral analytically, vorticity is taken to be linear as in (3). In addition, the integral is
decomposed according to the divergence theorem and is solved over the four faces of the tetrahedra. This expression
becomes:
(5)
This expression is used for nearby tetrahedra. Contributions of tetrahedra at an intermediate distance are computed
by 1-point or 5-point Gaussian quadrature. Contributions of tetrahedra further away are computed by an accelerated
calculation (Greengaard and Rokhlin, 1985).
Computation of Derivatives:
Since the vorticity is assumed to vary linearly over the element, the first derivatives will be constant over a single
element and the second derivatives will be zero. A higher order method to compute the derivatives was desired. A second
order method that expresses first and second order spatial derivatives across scattered points utilizes a polynomial fit of
the local vorticity using a least squares solution for all the tetrahedra which intersect a given node. For more detail, the
the authoritative version for attribution.
reader is referred to Marshall and Grant (1995). This is accomplished by expressing a component of the vorticity about a
desired node as:
(6)
and determining the constants a-k by a least squares fit to the values of ω−ωo at points in the neighborhood. x, y and
z are referenced to the local node. For tetrahedral methods, the elements connected to a given node are known and
therefore the nodal points can easily be determined. Additional points can be found through a search of neighboring
elements. After a sufficient number of points are found to maintain accuracy, a second order fit of the vorticity is
computed with the values of the derivatives simply computed at a given nodal point (x=y=z=0):
OCR for page 836
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 836
(7)
Evolution of the Vorticity Field:
In inviscid flow, the velocity field transports vorticity in the same way as a material element. This type of flow is
thus very well suited to Lagrangian mesh formulations, where the mesh points are transported by the velocity field.
However, when viscosity is present, vorticity may be transported by means other than advection by the velocity field
namely, via diffusion. Rather than introducing new ‘empty' points into the mesh onto which vorticity may diffuse, we
transport the existing mesh points with the sum of a diffusion velocity and the usual flow velocity. Thus the mesh points
tend to move from regions of larger vorticity magnitude to regions of lesser magnitude, according to the diffusive
transport by viscosity.
The concept of diffusion velocity for scalars (i.e. vorticity magnitude which can be written as one of the dependent
variables) is readily developed. Begin with the vorticity equation:
(8)
This equation states that the material change in vorticity is a function of vortex line stretching and diffusion due to
viscosity. Diffusion velocity based on the magnitude of the vorticity is defined as:
(9)
where Ω is defined as the scalar magnitude of the vorticity vector. This expression is inserted into equation (8) and
after some algebra, the result is:
(10)
Here, is the unit vector tangent to the vorticity vector defined as:
The mesh points are now transported according to:
(11)
and the vorticity is evolved according to (10).
The Lagrangian points were advected according to an Adams-Bashforth method to maintain second order accuracy
in time. Since time step remains constant:
(12)
Baldwin-Lomax Turbulence Model:
Since spatial resolution constraints due to the computational cost of the calculation prohibit direct numerical
simulation of high Reynolds number flow past the bodies of interest, a Baldwin-Lomax turbulence model was
implemented into the code. This turbulence model is summarized by Wilcox (1993). It is an algebraic eddy viscosity
model that was developed for use in computations where boundary layer properties such as boundary layer thickness and
edge velocity are difficult to determine. This is the case for the present unsteady flow cases. Since the flow is unsteady
and the computational points are effectively randomly scattered, these quantities would be nearly impossible to determine.
The Baldwin-Lomax model contains an inner layer and outer layer eddy viscosity:
Inner Layer:
(13)
Outer Layer:
(15)
(14) (16)
(17)
the authoritative version for attribution.
The closure coefficients are:
In order to determine y+, the friction velocity must be determined. In the present computations, the friction velocity
is assumed to be related to the surface vorticity value so that:
In the above formulation, y is the normal distance to the wall, ymax is the value of y where Fmax is found and Udif is
the maximum value of U for boundary layers. For free shear layers, Udif is the difference between the maximum velocity
in the layer and the value of U at y=ymax.
In determining whether to use the inner layer or the outer layer eddy viscosity, Fmax and ymax must be properly
determined. To properly determine these values, the computational points must be directly above the body surface. For
structured grids, these points are easily determined but for randomly spaced points, it becomes more difficult. A
OCR for page 837
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 837
subroutine was written to organize the volume points according to their nearest body point. In this manner, all points in
the volume are associated with a single body point that it is most nearly above. These sets of points are then sorted to
increase in the normal direction and then Fmax and ymax determined. If y is greater than ymax, the outer layer eddy viscosity
equation is used and if y is less than ymax, the inner formulation is used. Equations 14 and 15 are then modified using the
total viscosity defined as:
(18)
Vorticity Boundary Condition:
After the surface source and vortex sheet strengths are computed it is necessary to transfer the vorticity in the
infinitely thin sheets into the volume. This is accomplished by adjusting the surface nodal vorticity values to satisfy the
no-slip boundary condition exclusively. To do this, surface vortex sheet values must be minimized thus requiring:
(24)
The left-hand side is the volume integral of the first vorticity layer at node m that must balance the vorticity in the
infinitely thin layer expressed in the right hand side term. The vorticity is desired at body point m that can satisfy the no-
slip boundary condition. An iterative scheme is then used to set the vorticity. The surface velocity boundary conditions
are computed and the vorticity sheet strengths determined. The vorticity is then determined by equation 24. The surface
velocity boundary conditions are then re-computed and vorticity sheet strengths computed again. This iterative procedure
continues until the magnitude of the vortex sheet strength is below 0.01. This typically only requires 2–3 iterations and
converges quite rapidly. The reason for this is that the layer of vorticity is very thin so there is little difference in the
velocity generated by the infinitely thin sheet and the thin vorticity elements connected to the surface.
Initial Volume Vorticity Distribution, Euler Layer and Point Creation:
After the surface panels are defined, it is required to initialize the volume vorticity on a set of points. Nodal vorticity
values located at the surface body points are used to define the geometry. Additional points are then placed in layers
staggered over the body nodes and the surface panel centroids. To ensure that sufficient resolution of the boundary layer
vorticity is maintained close to the surface, a thin layer of “Euler layer” of fixed points is used. Typically, the field points
are located normal to the body nodal points in successive layers. This also allows for the computation of derivatives using
standard finite difference formulas as an alternative to the least squares approach summarized in equations 6 and 7.
Vorticity evolution for points in the Euler layer is performed using the Lagrangian form of 8–10 then interpolating the
values back on to the original point positions. On the first time step, the Euler Layer consists of seven sub-layers and ten
additional sub-layers of points are Lagrangian (advected by the local flow). Figure 2a shows a close-up view of the
typical initial distribution of points and subsequent cross-section of the tetrahedral mesh for a cone. Figure 2-b shows the
developed point distribution and mesh cross-section at a later time.
Figure 2: a) initial point distribution and b) point distribution at t=3.0.
On the first time step, an initial boundary layer thickness is assumed, which is 10% of a fully developed, turbulent
boundary layer on a flat plate. The boundary layer is assumed to be attached and of uniform thickness modeling an
impulsive start. The purpose of this is to more easily initialize the vorticity field with the vorticity remaining close to the
wall. The vorticity on the surface is determined from the boundary conditions and the vorticity in the volume is assumed
to exponentially decay. The vorticity values on the surface and in the boundary layer are then adjusted using an iterative
scheme to satisfy the no-slip and no-flux boundary conditions.
New Lagrangian points are continually created above the Euler layer. Extended panels are formed that are effectively
the authoritative version for attribution.
the outer faces of the tetrahedra formed in the Euler layer. Since the points in the Euler layer are not advected, the
tetrahedra in this layer can be explicitly set for all times. These extended panels on the outer surface of the Euler layer
have the same characteristics as the surface panel directly underneath. New Lagrangian points are formed by first
determining the closest point to the extended panel that is above the panel. To determine if the point is above the panel,
the centroid of the surface panel is determined and a radius constructed to include the three panel vertices. If a point lies
within the radius, its normal component is computed. The closest normal point within the radius is then found. If this
point is greater than a prescribed value, a new Lagrangian computational point is created. The new point is located in the
tetrahedron formed by the closest point and the three points connected to the nodes of the extended panel. The vorticity is
OCR for page 838
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 838
set to be the average of the four nodal vorticity values. If no Lagrangian point is found at all above the panel, a new point
is created above the panel at a prescribed distance away from the Euler layer and the vorticity is set by extrapolating the
vorticity from points within the Euler layer.
Figure 4: Axial and azimuthal vorticity component at
t=0.0, 1.0 and 2.0 for a columnar vortex with initial
Figure 3: rms error in the velocity calculation for various Gaussian distribution of vorticity. Symbols display the
accelerated calculation error bounds. exact solution.
Numerical Accuracy:
In order to define the numerical accuracy and grid convergence of the present method, some tests were performed.
First, the derivative computation is between first and second order accurate for irregularly spaced points. Test
computations were examined to compare the accuracy of the velocity calculation for increasing number of points, a test
using Hill's spherical vortex was conducted. This test has the advantage that an analytical solution of the velocity field is
known so direct comparison between the discretized solution and analytical solution can readily be made.
Figure 5: UUV and cone geometries
The computational points were randomly distributed inside a unit sphere with number of points between 1,000 and
32,000 used. In the accelerated calculation, different errors, ε, between 10−2 and 10−4 were examined. Figure 3 shows a
plot of the root mean square (rms) error as a function of the number of points for different error bounds. As a reference,
the vortex blob solution is plotted. As can be seen, the error for the vortex blobs are significantly higher compared with
the tetrahedral solution. As N is increased, the rms error for ε= 10−3 appears to decrease asymptotically to 10−2. For ε=10
−4, the rms error is only slightly improved.
Validation computations were also performed to verify the vorticity diffusion algorithm by comparing the analytical
solution for a columnar vortex with the present solution. Figure 4 shows the axial and azimuthal vorticity components at
t=0.0, 1.0 and 2.0 for the computed solution and the exact solution. As can be seen, the computed solution agrees very
will with the exact one.
UUV and Cone Definition:
For all test problems, non-dimensional geometric quantities were used. Maximum freestream velocity is defined as
1.0. UUV body length is also 1.0. A time of 1.0 is defined as the time it takes the freestream to travel one body length.
the authoritative version for attribution.
Reynolds number can then be set as an independent variable and for these cases was set to 1,000,000 for flow past a UUV
and 50,000 for flow past a cone.
The UUV was defined using 550 surface points resulting in 1100 panels (Figure 5a). The cone was defined using
442 points resulting in the formation of 880 panels (see Figure 5b). A 45° cone angle was used with diameter equal to
two, length equal to one and flat stern face. The cone nose was placed at the origin of the global Cartesian coordinate
system and the nose of the docking UUV was placed at distances (XUUV) of 0.5, 0.25 and 0.125 relative to the origin.
For all test cases, the flow is impulsively started from a freestream velocity of 0.0 to 1.0. During this initial startup,
OCR for page 839
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 839
the vorticity is diffused onto the initial layers of points forming a fully attached boundary layer over the cone as well as
the UUV as was described earlier. As time progresses and the flow develops, a separated region forms behind the cone
characterized by a cohesive vortex ring structure. This vortex forms close to the surface of the cone and the resulting
wake structure impacts the downstream UUV.
Figure 6: Instantaneous tetrahedral grid cross-section in
the x-y plane at t=2.0.
Figure 8: Turbulent boundary layer profiles comparing
the current computations with empirical flat plate results
from Spalding and Cole.
Figure 7: Instantaneous vector plots in the x-y plane of Figure 9: Same as Figure 8 except scales adjusted to
the turbulent flow past a UUV. The vectors are colored highlight near wall region.
based on the vorticity magnitude.
Test runs were conducted at both the Waterways Experimental Station Cray C-90 and at the NAVOCEANO Cray
C-90. Average run times required approximately 36 CPU hours.
RESULTS:
A cross-section of the tetrahedral mesh depicting the instantaneous flow past a UUV at zero angle of attack for a
Reynolds number of 1,000,000 is shown in Figure 6. This figure illustrates the resolution of the grid points used to
compute the unsteady turbulent flow. Figure 7 shows snap
Figure 10: Turbulent intensities in velocity components (u', v' and w') compared with flat plate empirical results from
the authoritative version for attribution.
Klebanoff.
OCR for page 840
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 840
shots of the vector plots taken at x/l=0.7 for t=2.0, 2.1 and 3.0. The vectors are colored based on the vorticity magnitude
with the scale shown to the right. Maximum vorticity magnitudes of 100.0 were chosen to illustrate the boundary layer
flow. Maximum surface vorticity values exceeded 2,000. These plots appear very similar. Upon closer inspection,
however, there are qualitatively small changes in the overall flow as the unsteady vorticity field is advected over the
UUV. These plots illustrate the unsteady attached flow over the UUV for high Reynolds numbers.
Boundary layer axial velocities over the UUV were time averaged between t=2.0 and t=4.0 to obtain mean velocity
profiles and the turbulent intensities in the velocity components (u', v' and w'). Mean axial velocity profiles for
Re=1,000,000 are plotted in Figure 8 for the current computations and for empirical results obtained from flat plate
boundary layers for Spalding and Law of the wake corrections of Coles (see White, 10). As can be seen, the computed
velocity profile agrees very well with empirical solutions from Spalding and Coles. Even in the inner layer, the agreement
is good except for slightly lower velocities at y/δ between 0.01 and 0.04.
Figure 9 shows the same mean velocity profile with the scales compressed to highlight the near wall behavior. At the
wall, the computations agree very well. From y/δ=0.01−0.04, it appears that the computational solution bows out
somewhat. Good agreement is seen thereafter from y/δ=0.04−0.1.
Figure 10 shows plots of the turbulent intensities (u', v' and w') as a function of distance with comparisons with these
statistics obtained by Klebanoff (see Hinze, 1975). Since a turbulence model is used, the turbulent intensity statistics were
expected to be minimal. The computations, however, are showing fluctuations in all three velocity components due to the
inherent unsteady computations of the present method. u' appears to show the best agreement although it does not reach
the same maximum near the wall as Klebanoff's data. v' and w' show poorer agreement. This is likely due to the lack of
spatial resolution of the surface. The fact that the turbulent intensities are non-zero and the fact that u' follows the same
trends as the experimental results demonstrates the potential of the current code to represent turbulent flow without a
turbulence model. Of course, the spatial resolution would need to be greatly increased to properly accomplish this.
Unsteady flow past a 45° cone was also computed and results of the unsteady flow field can be seen in Figure 11.
Here, vector plots, representing the instantaneous velocity field, are colored based on the z-vorticity component with the x-
y plane displayed. The surface is shaded based on the vorticity magnitude. This plot compares the instantaneous unsteady
flow at t=2.0, 3.0 and 4.0 for Re=1,000 and Re= 50,000. For Re=1,000, the flow was assumed laminar and for
Re=50,000, the turbulence model was used. In both cases, a ring vortex is produced in the wake of the cone and can be
seen to grow in size as time is increased. The laminar flow case shows a much more coherent, well-defined vortex
compared with the turbulent flow case.
Quantitative comparisons of this flow condition are shown in Figure 13 for laminar and turbulent flow computations.
The z-vorticity component in the wake was computed at one cone diameter downstream of the stern of the cone
(corresponding to the location of the vortex core). The computational data shows the instantaneous velocity field at t=4.0.
As can be seen, data for Re=1,000 show significantly higher vorticity with peak values on the order of 20. Turbulent flow
vorticity values for Re=50,000 are much more modest by comparison with values approaching 10. In addition, the
vorticity for the laminar flow case appears to be
Figure 12: Z-Vorticity distribution at t=4.0 in the wake
of the cone at x=2.0 (relative to the cone apex) for
Re=1,000 (laminar flow) and Re=50,000 (turbulent flow).
the authoritative version for attribution.
Figure 11: Unsteady flow development past a cone for
the laminar flow case (Re=1,000) and the turbulent flow
case (Re=50,000). Velocity vectors are colored based on
the wake z-vorticity and the surfaces are colored based
on the vorticity magnitude.
OCR for page 841
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 841
more concentrated within a well-defined core. The turbulent flow case appears to be more spread out.
Figure 13: Wake velocity distributions for a 45° cone 0.88 diameters downstream of the cone stern for Re=50,000 and
Re=1,000 at t=5.0. The mean wake profile from Calvert is shown as a reference.
Streamwise velocity component data are shown in Figure 12 for laminar and turbulent flow computations. The
streamwise velocities in the wake were computed at 0.88 cone diameters downstream of the stern of the cone. This
corresponds to the same location as Calvert (1967) for his experimental wake velocity data for flow past a cone. The
computational data shows the instantaneous velocity field at t=5.0. Data from Calvert are time averaged. As can be seen,
data for Re=1,000 show significantly stronger reverse flow velocities along the centerline of the cone with reverse flows
approaching −1.5. Also, maximum flow velocities are on the order of 1.5 near the shear layer where the freestream meets
the vortex wake. Turbulent flow computations produced milder flow velocities with maximum streamwise velocities of
1.1 seen in the shear layer and reverse flow velocities reaching −1.0. These data also compare well with experimental
results. The initial decrease in streamwise velocity (y/d=0.3–0.6) compares quite favorably as do maximum streamwise
velocities. Experimental results show reverse flow velocities on the order of −0.5. Since we are comparing instantaneous
computational results with time averaged, however, the computational results appear highly encouraging.
UUV Docking Simulation:
The docking cone was defined using 442 points resulting in the formation of 880 panels (see Figure 14). A 45° cone
angle was used with 0.0954 diameter and the stern face of the cone was flat. The nose of the cone was placed at the origin
of the global Cartesian coordinate system and the nose of the docking UUV was placed at distances (XUUV) of 0.5, 0.25
and 0.125 relative to the origin.
Figure 15 shows the initial flow development for the test case of the UUV at a position 0.125 downstream of the
origin (slightly less than one cone diameter downstream of the cone stern). The surface displays the source/vortex panels
with the surface colored based on the y component of surface vorticity (from ±500). The vectors in the field represent the
velocity at the tetrahedra centroids and are colored based on the streamwise component of the velocity with red vectors
displaying maximum velocities over two times freestream values and blue vectors displaying reverse flow velocities. At
t=0.0, the effective attached boundary layer flow can be seen. Notice the high velocities at the top and bottom of the cone.
Immediately aft of the cone, the velocity vectors point toward the axis of the cone with zero streamwise velocity
component. The initially structured mesh points over the UUV display a well defined boundary layer produced by
diffusing the initial vorticity onto the points over the surface. After 10 time steps (t=0.05), the separated vortex ring
structure can clearly be seen aft of the cone. By t=0.1, this vortex impacts the leading edge of the UUV. Notice a slight
reduction in surface vorticity by the decrease in redness near the leading edge. By t=0.15, the vortex in the wake becomes
larger and the surface vorticity becomes redder indicating increased vorticity. By t= 0.2, impact of the vortex ring
structure on the nose of the UUV
the authoritative version for attribution.
Figure 15: Velocity vector plots at the tetrahedra
centroids depicting the unsteady flow development past a
Figure 14: UUV, Cone and tandem UUV/Cone geometry. cone/docking UUV from t=0.0–0.5. Re=10,000,
XUUV=0.125. Vectors are colored based on the
streamwise velocity component and the surface is
colored based on the y vorticity components.
OCR for page 842
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 842
creates an asymmetry. In addition, notice the blue surface coloration indicating opposite signed vorticity production at the
nose of the UUV. Beginning at t=0.25, periodic asymmetry can be seen in the vortex structure between the cone and the
UUV. Finally, notice the increase in boundary layer thickness and the formation of smaller vortex structures as the
vorticity from the cohesive vortex ring is ingested into the UUV boundary layer.
Figure 17: Velocity vector plot comparison of flow
Figure 16: Velocity vector plot comparison of flow
along the centerline in the x-y plane for XUUV=0.2 and
along the centerline in the x-y plane for XUUV=0.125,
the UUV center vertically displaced—0.05. The surfaces
0.25 and 0.5 at t =1.0.
are colored based on the surface pressure and the velocity
vectors are colored based on the z-vorticity component.
Figures 16 shows comparisons of the flow for XUUV= 0.125, 0.25 and 0.5 at time 1.0. For XUUV=0.125, the flow
from top to bottom is mostly symmetric although the alternate variation in vortex size does appear to result in increased
flow over the upper surface by t=1.5. For XUUV=0.25 and 0.5, the UUV is sufficiently downstream of the cone so a
definite oscillation due to changes in wake vorticity shedding occurs. This wake oscillation results in different impacts on
the trailing UUV. This results in an overall thickening of the UUV boundary layer with areas of vorticity accumulations,
which include the vorticity produced by the cone.
In Figure 17, the UUV is displaced a vertical distance of 0.05 below the axial centerline of the cone. The cone is at
0° pitch angle and the XUUV=0.2. Again, the surface is colored based on the surface pressure and the velocity vectors are
colored on the z-vorticity component. In this case, the vortex structure on the upper portion of the cone is allowed to
freely advect over the UUV while the vortex on the lower portion of the cone impacts the UUV. At t=0.1, the cone vortex
wake appears fairly symmetric. By t=0.2, the initial portion of the wake impacts the UUV. This results in the formation of
a vortex structure near the leading edge. By t= 0.3, this wake induced vortex begins to advect over the UUV. A low-
pressure region is indicated immediately downstream of this vortex. By t=0.4 and 0.5, it appears that the entire wake
produced by the cone advects over the upper surface of the UUV and combines with the boundary layer vorticity of the
UUV. There did not appear to be any portion of the cone wake that advected over the lower portion of the UUV.
Figure 18: Surface pressure on the UUV as a function of
Figure 19: Same as Figure 8 except XUUV=0.5.
time for z=0.0, upper surface and x=0.0, 0.06 and 0.5
corresponding to nose stagnation point, point of
maximum UUV radius and UUV body center
respectively. XUUV= 0.125.
Figure 18 shows the point pressure distributions as a function of time for the XUUV=0.125 test case. The point
locations correspond to UUV nose stagnation point (x=y=z=0.0), point of maximum radius on the UUV (x=0.06,
y=0.045, z=0.0) and the center of the UUV (x=0.5, y=0.045, z=0.0). This plot demonstrates the transient nature of the
surface pressure and highlights the impact of the cone wake on the UUV. The pressure at the nose shows an initial
stagnation pressure of 1.0. By t=0.1, the pressure decreases to 0.0. By this time, the vortex ring structure is well formed in
the wake. At t=0.25, there is an increase in pressure which correlates with full impact of the vortex on the UUV nose. As
the authoritative version for attribution.
time progresses, there are fluctuations in the pressure.
Pressure at maximum radius (x=0.06) shows a different behavior. Initially, Cp values of −0.3 are seen before impact
of the cone wake. There appears to be an increase in
OCR for page 843
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 843
Cp to 0.0 at t=0.35 which correlates with boundary layer ingestion of the wake vorticity. Additional fluctuations correlate
with the ingestion of accumulated vorticity in the boundary layer. At the center of the UUV (x=0.5), Cp values remain at
approximately 0.0 until t=1.125 where a sharp decrease in pressure is observed. This Cp spike is quite transient and
correlates with a locally strong vortex in the boundary layer.
Figure 20: Unsteady force and moment coefficients for Figure 21: Unsteady force and moment coefficients for
the UUV with the origin at 0.125 the UUV with the origin at 0.5
Figure 19 shows the point pressure distributions as a function of time for the XUUV=0.5 test case. Nose pressure
shows stagnation Cp values of 1.0 to t=0.45. This is followed by a sharp decrease in Cp to 0.0 by t=0.5 followed by an
immediate increase back to Cp=0.8. Cp values slowly decrease with fluctuations in pressure seen. Pressure at the
maximum radius shows steady Cp values of −0.5 up to the point of wake impact on the UUV at t=0.5. As can be seen, at
t=0.55 a sharp increase in Cp to 0.2 is seen followed by a decrease in Cp back to −0.5 by t=0.75. Pressure then remains
approximately constant until just after t=1.0 where additional fluctuations are seen as the wake becomes ingested in the
UUV boundary layer. At the UUV center, Cp values of 0.0 persist out to t=0.7 where a decrease to −0.3 is seen at t= 0.75.
This correlates with the passage of a vortex. As time is increased and the vorticity is advected downstream, Cp values
recover back to 0.0 with fluctuations in Cp seen thereafter.
Figures 20 and 21 show the UUV forces and moments as a function of time for UUV positions of 0.125 and 0.5
respectively. In these cases, the forces are nondimensionalized by the product of the freestream dynamic pressure and
wetted surface area; the moments are nondimensionalized by the product of the freestream dynamic pressure, wetted
surface area and vehicle length with the moment taken about the vehicle center (x=0.5, y=0.0 and z =0.0. For all three
cases, the z component of the force is negligible due to the assumption of symmetry about the y-axis. As a consequence,
the y-moment component is negligible as well as the rolling moment about the x-axis.
For XUUV=0.125, drag coefficient appears to remain level at 0.025 out to t=1.3. After that, there is a rise and
fluctuation to approximately 0.04 followed by another decrease. Normal force (Cfy) fluctuates about the zero value until
just before t=1.0. Afterward, there are several fluctuations with a large transient peak in Cfy of 0.3 at t=1.35. This is
followed by another peak of 0.1 at t=1.5 before the forces diminish again. The pitching moment (Cmz) shows small
values until t=0.75 where a decrease in moment to −0.01 is seen. Afterward, there appears a definite fluctuation in the
UUV moment with peak values on the order of 0.02 but are very transient in nature. It is interesting to note that even
though the wake impacts the UUV by t=0.15, force fluctuations are not seen until t=0.75 at the earliest.
At the furthest downstream location of 0.5, it appears that all fluctuations in the UUV forces and moments are
diminished. The only fluctuation appears at t=0.5 just as the cone wake impacts the UUV. After the cone vorticity is
ingested in the boundary layer, vehicle normal force and pitching moment return to effectively zero values and the drag
coefficient remains constant at 0.025.
Figure 22 shows the unsteady forces and moments on the UUV for the cases for XUUV=0.2 and YUUV=−0.1. Drag
force remains fairly constant at 0.02 and side force remains negligible. Notice that the normal force becomes negative by t
the authoritative version for attribution.
OCR for page 844
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 844
=0.1 with values remaining fairly constant up to t=1.2 at Cfz=0.01. After that time a decrease in normal force to Cfy=−0.1
is seen out to the end of the run. Pitching moment shows the effect of the displaced wake as well. Initially, pitch moment
becomes positive by t=0.4 with Cmz values on the order of 0.01. The increase in normal force correlates with an increase
in pitch moment suggesting a nose down pitch moment with downward forces proximal to the UUV nose.
Figure 23: Bode plots displaying the gains in vertical
position due to pitch moment and vertical force.
Figure 22: Unsteady force and moment coefficients for
the UUV with the origin at (0.25, −0.10, 0.0).
DISCUSSION:
Unsteady Wake Development:
The velocity vector plots displayed the initial unsteady flow development aft of the docking cone and the resultant
impact of the wake on the docking UUV. After the flow travels approximately one cone length, the initial vortex ring
structure aft of the cone is formed. This ring structure then becomes elongated and appears to form an elliptically shaped
vortex ring by t=0.2. As long as the UUV is not in close proximity (e.g. XUUV=0.5), the major axis of the vortex ring
appears to reach a maximum size about twice that of the minor axis. For the UUV in close proximity, the cross-section of
the vortex ring remains circular with fluctuations in the local diameter of the ring vortex as the flow develops.
The presence of the UUV appeared to alter the frequency of the force fluctuations on the UUV. For XUUV= 0.125
and 0.25, high frequency fluctuations in drag and cone pitching moment were seen. For the far downstream UUV case,
however, a relatively well-defined Strouhal frequency of 0.225 was observed. This suggests that the UUV alters the
vorticity shedding process of the cone as it comes in close proximity. As the flow continues to progress, asymmetry in the
vortex ring structure results in the shedding of ‘chunks' of vorticity. These vorticity ‘chunks' do not appear as well formed
vortex rings, rather they may more appropriately be described as convoluted hairpin vortices. It is this vorticity that
impacts the UUV and alters the local pressure distribution. For the far downstream cases, there appears to be little
influence of the wake on the integrated forces and moments although significant variations in surface pressure were
observed. As the UUV gets closer to the docking cone, however, significant unsteady pressure and loads are produced.
There is no longer a stagnation point at the nose of the UUV and the low pressure at maximum UUV radius normally seen
is diminished. The impacting vorticity leaves a definite signature in the pressure distribution and generates significant
normal force and pitching moment.
Full-scale UUV Recovery:
A full-scale UUV is on the order of 6.096 m (240”) long with 53.34 cm (21”) diameter. The cone diameter is of the
order 55.88 cm (22”). It is envisioned that the UUV will dock with the cone for a submarine forward velocity of 2.572 m/
sec (5 knots). Using these numbers, a force coefficient of 0.01 corresponds to a force of 280 N (63 lbs.) and a moment
coefficient of 0.01 corresponds to 1625 N-m (1200 ft-lbs.). A non-dimensional time unit of 1.0 corresponds to a full-scale
time of 2.11 seconds. Fluctuations with a non-dimensional period of 1.0 therefore correspond to a full-scale frequency of
0.474 Hz.
the authoritative version for attribution.
Figure 23 shows Bode plots for the 21” UUV displaying the gain in vertical position for inputs of pitching moment
and vertical force. These plots effectively demonstrate the response of the vehicle to unsteady loads. It
OCR for page 845
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 845
assumes that a sinusoidal forcing function is input with given amplitude. Actual displacement (sinusoidal) can be related
to the gain by the equation:
(25)
where y is the displacement, G is the gain and F is the forcing function (either force or pitching moment). These
plots show some revealing characteristics. First of all, there is a peak frequency at 0.02 Hz. This suggests that the vehicle
will be unstable to forcing functions with a period of 50 seconds. As frequency is increased, the gains decrease
significantly so that by 1 Hz, the gain is only on the order of 10−5 for the pitching moment. At even higher frequencies,
the gain drops off sharply so that, in effect, the vehicle doesn't even respond to the forces.
As could be seen from Figures 9–11, the high frequency vortex shedding resulted in force fluctuation periods on the
order of 0.1 seconds for resulting frequencies of 10 Hz. The corresponding gain is on the order of 10−8 resulting in
negligible change in UUV trajectory. Strouhal shedding frequencies on the order of 0.75 Hz were observed from the
moment fluctuations in Figure 11. The corresponding gain in pitching moment is 1.356x10−5 and is 1.334x10−3 in vertical
force. Assuming a forcing function that produces 1625 N-m (1200 ft-lbs.) of pitching moment and 280 N (63 lbs) of
vertical force, the resulting change in trajectory is on the order of 9.1 mm (0.36 inches). This change should be corrected
during a docking procedure and is easily within the window of error.
Although the vortex shedding off of the cone will likely have minimal impact on the trajectory of the UUV, motion
of the submarine or the presence of ocean currents will cause the cone wake to vary as well. A sinusoidal wave pattern
could be established on the wake of the cone. For moderate ocean currents, and temporal variation with 10 second period
superimposed on the cone wake would not be unreasonable. This wave pattern could result in pitch and moment
fluctuations on the order of 0.1 Hz. The bode plots shows a gain in pitch moment of 3.8x10−3 and gain in vertical force of
8.58x10−2. For a 1625 N-m (1200 ft-lbs.) pitch moment amplitude and 280 N (63 lb) vertical force amplitude, this would
result in changes in vehicle trajectory on the order of 1.37 meters (4.5 feet). This altered trajectory could be significant in
terms of docking the UUV. Future flow calculations should focus on these types of phenomena.
CONCLUSIONS:
A novel Lagrangian vorticity method has been presented to compute the unsteady hydrodynamics associated with
naval unmanned undersea vehicles. Since the flow is incompressible, the pressure term is not required for solution leaving
a velocity-vorticity formulation. In fact, since the velocity is a direct integral quantity of the vorticity, this method
demonstrates that the unsteady flow can be described by the vorticity alone. The Lagrangian nature of the calculation
requires only that the diffusion term be solved explicitly. The advection term is automatically included since the points
are moved with the local flow. The diffusion velocity concept is also used to move the points into regions of zero
vorticity. This avoids the difficulty of establishing empty points to diffuse the vorticity onto. The diffusion equation was
modified accordingly and solved at each time step. Effects of grid resolution and validation of the velocity calculation and
diffusion algorithm were conducted by comparing computational results for analytical solutions for Hill's spherical vortex
and a columnar vortex.
The unsteady turbulent flow past a UUV at Re= 1,000,000 and the unsteady flow development in the wake of the
cone were both investigated to demonstrate the effectiveness of the current method. Mean turbulent boundary layer
velocity profiles agreed quite well with empirical results from Spalding and Cole. Turbulent fluctuations agreed
surprisingly well for u' but not as well for v' and w'. This was expected due to lack of surface resolution. The fact that
turbulent quantities are computed at all demonstrates the potential of this method to compute turbulence directly.
Unsteady flow past a cone showed the vortex ring structure produced. Initially, the vortex ring was symmetric but
developed asymmetries as the flow developed. The laminar flow case showed a much more coherent vortex compared
with the turbulent vortex. This is exactly as intended and was the reason for implementing a turbulence model.
Consequently, the turbulent flow results for the wake velocity profiles were in much better agreement with experimental
test cases.
Simulation of UUV recovery hydrodynamics displayed some interesting flow field phenomena. The unsteady flow
development in the wake of the cone showed the vortex ring structure produced. Initially, the vortex ring was symmetric
but developed asymmetries as the flow developed. As the asymmetry was manifested, shedding of vorticity was observed.
The shed vorticity did not appear as well-formed vortex ring but as convoluted hairpin vortex structures. Proximity of the
UUV affected the way in which the vorticity was shed. For cases where the UUV was close to the UUV, the separated
vorticity region between the UUV and the cone maintained coherent vortex ring structures whose asymmetry varied with
time. As the UUV was placed downstream, the vortex ring structure became elongated and elliptically shaped. Here,
definite ‘chunks' of vorticity were shed which impacted the UUV.
As the wake of the cone impacted the UUV, significant transients in local pressure were observed. There was a
reduction in stagnation pressure at the nose of the UUV and increase in pressure at the point of maximum radius of the
the authoritative version for attribution.
UUV. Pressure fluctuations were observed as the vorticity dominated wake impacted the nose of the UUV. Calculations
at mid-body showed fluctuations in pressure as well. Even though there was large pressure fluctuations, integrated force
computations showed little fluctuations during initial wake impact. Only as the flow became developed were significant
force fluctuations seen. Also, the magnitude of the force fluctuations rapidly diminished, as the UUV was placed
downstream.
Vertical placement of the UUV produced altered wake structure and unsteady loading on the UUV. Placing the UUV
below the cone altered the vortex wake interaction with the UUV and also produced negative normal forces and nose
down pitch moment. For this case, sinusoidal fluctuations in
OCR for page 846
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 846
pitch moment on the order of 1200 ft-lbs and vertical forces of 60 lbs were observed. This was likely due to the overall
momentum deficit in the wake resulting in higher pressure on the UUV upper surface.
An analysis was performed to determine the resultant changes in UUV trajectory for a full-scale UUV. It was shown
that the fluctuations due to vortex shedding were sufficiently high to have minimal impact on the UUV trajectory. The
presence of the wake however and modifications of the wake due to submarine motion or ocean currents may adversely
impact the UUV trajectory. A wave-like wake structure with 10-second period was shown to alter the trajectory of the
UUV by 1.37 meters (4.5 feet). Future calculations should focus on the alteration of the docking cone wake due to
submarine motion and its resulting impact on the UUV.
Now that a turbulence model has been introduced into the code, it is possible to treat a wide range of engineering
problems of naval interest. Current work includes examination of unsteady flows due to bow thrusters for low speed
control, submarine maneuvering problems and UUV recovery problems. Ongoing work will be presented in future papers.
ACKNOWLEDGMENTS:
This work was sponsored by the United Kingdom Defence Evaluation and Research Agency under FMS case
numbers UK-P-GVK and UK-P-GUW, Richard Breward, DERA Program Manager, Simon Corfield, Technical Leader
and Mark Harbige Project Engineer and by the Office of Naval Research under contract number 99WX20012, Dr. Spiro
Lekoudis Program Manager.
REFERENCES:
Borouchaki, H., Lo, S.H., “Fast Delaunay Triangulation in Three Dimensions,” Computer Methods in Applied Mechanics and Engineering, Vol 128
(1995), pp. 153–167.
Calvert, J.R., “Experiments on the low-speed flow past cones,” Journal of Fluid Mechanics, Vol. 27, part 2, 1967, pp. 273–289.
Chorin, A.J., “Numerical Study of Slightly Viscous Flow,” Journal of Fluid Mechanics, V. 57, 1973, pp. 785–796.
Gharakhani, A., Ghoniem, A.F., “Simulation of Three-Dimensional Internal Flows by the Random Vortex and Boundary Element Methods,” European
Series in Applied and Industrial Mathematics, Vortex Flows and Related Numerical Methods II, ed. Gagnon, Y., Cottet, G.-H, Dritschel, D.G.,
Ghoniem, A.F., Meiburg, E., 1996 (http://www.emath.fr/proc/Vol.1/).
Greengard, L, Rokhlin, V., “A Fast Algorithm for Particle Simulations,” Mathematics of Computations, Vol. 47, pp. 387–398, 1987.
Hinze, J.O. Turbulence, 2nd Edition, McGraw-Hill, 1975.
Hoerner, S.F., Fluid-Dynamic Drag, Hoerner Fluid Dynamics, Bricktown, NJ, 1965.
Huyer. S.A., Grant, J.R., “Solution of the Two-Dimensional Vorticity Equation on a Lagrangian Mesh,” AIAA Journal, Vol. 38, No. 5, May 2000, pp.
774–783.
Koumoutsakos, P., Leonard, A. “High-Resolution Simulations of the Flow Around an Impulsively Started Cylinder Using Vortex Methods”, Journal of
Fluid Mechanics, Vol. 296, 1995, pp. 1–38.
Marshall, J.S., Grant, J.R., “A Lagrangian Collocation Method for Vorticity Transport in Viscous Fluid Flows,” proceedings from Forum on Vortex
Methods for Engineering Applications, Albuquerque, NM, Feb 22–24 1995.
Russo, G., Strain, J.A., “Fast Triangulated Vortex Methods for the 2D Euler Equations,” Journal of Computational Physics, Vol 111, pp 291–323, 1994.
Sethian, J.A., Ghoniem, A.F., “Validation study of vortex methods.” Journal of Computational Physics, Vol. 74, 1988, pp. 283–317.
Strickland, J.H., Baty, R.S., “An overview of fast multipole methods.” Sandia National Laboratories Report SAND95–2405, Albuquerque, NM,
November, 1995.
Subramaniam, S. “A New Mesh-Free Vortex Method”, Ph. D. dissertation, Dept. of Mechanical Engineering, Florida State University, 1996.
Uhlman, J.S., Grant, J.R., “A New Method for the Implementation of Boundary Conditions in the Discrete Vortex Element Method,” ASME 1993
Fluids Engineering Spring Meeting, Washington, D.C., June 1993.
White, F.M, Viscous Fluid Flow, 1st Edition, McGraw-Hill, 1974
Wilcox, D.C. Turbulence Modeling for CFD, Griffin Printing, Glendale, CA., 1993.
Zienkiewicz, O.C., The Finite Element Method, 3rd edition, McGrawHill Publishing Co., 1977, pp 165–169.
the authoritative version for attribution.
OCR for page 847
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
SIMULATION OF UUV RECOVERY HYDRODYNAMICS 847
DISCUSSION
U.Bulgarelli
Instituto Nazionale per Studi ed
Esperienze di Architettura Navale, Italy
In the equation of vorticity, how do you compute the coefficient u at t=0 for starting the computation?
AUTHOR'S REPLY
At t=0, an impulsive start is assumed. In this case, the freestream velocity is accelerated to 1.0 in one time step. The
vorticity is distributed on an initial set of points in the boundary layer with the point distribution normal to the surface.
The thickness of this initial boundary layer is assumed to be 10% of the final thickness of a flat plate boundary layer. On
the first time step, the no slip, no-flux boundary conditions are satisfied with the surface vortex sheet. The surface
vorticity is then set so that ωdV=γdA where ω is the surface vorticity, γ is the surface vortex sheet strength, dV is the
volume of the elements connected to the surface node and dA is the area of the panels connected to the surface node. The
vorticity is then assumed to exponentially decay through the initial layer with a maximum vorticity at the surface. Since
the distribution of vorticity will affect the velocity at the surface, an iterative method was used to re-calculate the surface
vortex sheet, surface vorticity and layer vorticity. By the end of this iterative process, the vorticity and velocity of each
point are known.
the authoritative version for attribution.
Representative terms from entire chapter:
flow past
