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

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

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:

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):

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

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

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,

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.

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.

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.

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

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.

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

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

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.

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.