Twenty-Fourth Symposium on Naval Hydrodynamics (2003)

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.

24th Symposium on Naval Hydrodynamics Fukuoka, Japan, 8-13 July 2002 Unstructured Nonlinear Free Surface Simulations for the Fully-Appended DTMB Model 5415 Series Hull Including Rotating Propulsors Clarence 0. E. Burg, Kidambi Sreenivas, Daniel G. Hyams, Brent Mitchell (Computational Simulation and Design Center, Mississippi State, MS) ABSTRACT Nonlinear free surface simulations around realistic ge- ometries, such as the DTMB Model 5415 Series hull, are a necessary step to achieve the goal of simulation of maneuvering surface vessels. As part of this effort, re- searchers in the Computational Simulation and Design Center at Mississippi State University have developed and parallelized a three-dimensional unstructured code, U2NCLE, which solves the incompressible Reynolds- averaged Navier-Stokes equations using a node-based flux differencing finite volume method with Roe-averaged variables. The simulations produced by this code real- istically capture turbulent flow and vortices arising from bulbous bows and tips of propulsors and rudders and can model the rotation of the propulsors accurately. Because the code uses unstructured grids, the code is applicable to a wide range of geometries, including those with shafts and struts, and can solve the flow through small gaps such as those between the rudder and the boot. Current re- search efforts include the addition of nonlinear free sur- face capabilities to this code. A nonlinear free surface is obtained by solving the kinematic free surface boundary condition, imposing the resulting hydrostatic pressure dis- tribution onto the boundary for the Navier-Stokes solver and moving the grid to match the free surface while con- forming to any geometry that intersects the free surface. The current free surface implementation is loosely cou- pled with the Navier-Stokes solver and makes the assump- tion that the flow is steady-state. INTRODUCTION Modeling the free surface generated by surface ships and by submerged vessels near the water surface is important for proper understanding of the flow around these ves- sels. In particular, the free surface changes the resistance of the vessels through the water by changing the pres- sure distribution and the wetted surface area. The free surface affects the location and magnitude of the vortices that originate from various locations, including the bow, appendages and propulsors, which can greatly affect the performance of the propulsors. The wave signatures of ships and submarines is also of importance. Surface ves- sels under investigation by ship designers attempt to miti- gate the negative influences of the free surface interaction on cavitation, power requirements and wave signature and attempt to use the free surface to improve the performance in other ways. Because of the novelty of many of these designs, which often are not direct extensions of current designs, extensive sets of experimental data may not yet exist, and designers are more apt to rely on data from nu- merical simulations. Furthermore, the effects on the ship's performance caused by maneuvering through a change in direction, or through various sea states, especially in re- gards to sea-keeping and cavitation in littoral water are of importance in the designs and can not be readily tested by experiment. A long-term goal of the research into numerical simulations is to develop the ability to study the perfor- mance of a full ship design, including the interactions of the various appendages, sonar domes, rudders, shafts and propulsors, as the ship maneuvers in response to changes in the settings of the rudders and propulsors. To accom- plish this goal, efficient and accurate free surface simula- tions are needed. Other necessary components to achieve this goal include capabilities to monitor the onset of incip- ient cavitation, to perform simulation of several vessels in motion relative to one another and to adapt the grid to cap- ture vortices in the water. Currently, the implicit unstructured code being developed by researchers in the Computational Simu- lation and Design Center at Mississippi State Univer- sity solves the three-dimensional, incompressible Navier- Stokes equations, via an edge-based, flux-differencing fi- nite volume method with Roe-averaged variables. Hyams (2000a, 2000b) successfully parallelized the code and has ported the code to a wide variety of high-performance ma- chines. The unstructured grid generator (Marcum, 1998) has the capability of building highly stretched high-aspect ratio grids near viscous surfaces, which include higher or- der elements such as pyramids and prisms. Other areas of code development include realistic rudder and propulsor- induced motion and full-scale Reynolds number simula- tion.

At each time step, the nonlinear free surface al- modeled under these assumptions. gorithm solves the kinematic free surface equation BY BY ot + (u—V=) ~ —(v - Vy) + (w—Vz) ~ = 0 (1) where Y = Y(x,z,t) is the free surface defined as a single-valued function over the xz plane, (u,v,w) are the velocity components in the coordinate directions and (V~, Vy7 Vz) are the grid velocities. A flow-throughbound- ary condition for the free surface based on characteristic variable boundary conditions is used, with the pressure on the free surface set to P = ~ where P is the pressure and Fr is the Froude number. After several time steps, the grid is moved to match the free surface while conforming to the geometry, using a three dimensional extension of Farhat's torsional spring method (Farhat, 19981. This grid movement algorithm is quite robust, allowing for moder- ate to extreme distortions as required by the free surface simulations. As the solution converges, the flow along the free surface becomes tangent to the free surface, and the nonlinear free surface solution is obtained. Much of this work is an extrapolation of the free surface algorithm used within the three-dimensional struc- tured code UNCLE also developed at Mississippi State University. Beddhu (1994)developed the structured free surface solver, using a modified artificial compressibility formulation. This algorithm was applied to steady and unsteady flow around the Wigley hull (Beddhu, 1998a), to the barehull model 5415 series hull (Beddhu, 1998b), to the Series 60 CB = 0.6 ship (Beddhu, 1998c) and to a more detailed stern analysis for the model 5415 series hull (Beddhu, 1999), and these results were presented at the Gothenburg conference (Beddhu, 2000~. Initial veri- fication and validation exercises of the unstructured non- linear free surface algorithm has been presented by Burg (2002), which includes a grid refinement study for flow over a submerged NACA0012 hydrofoil, flow around in- viscid and viscous Wigley hulls, and flow around the bare- hull model 5415 series hull. In this paper, results for the DTMB Model 5415 Series hull, which has a transom stern, are presented. The numerical results are compared against experimental re- sults which include the hull profile for the barehull 5415 and the stern region topologies for the unpowered and powered fully-appended 5415. Because of the use of un- structured grids, the actual geometry with the shafts, struts, rudder, the gap between the rudder and the boot and the DTMB propellers 4876 and 4877 is used to generate the grid. For the unpowered case, the propulsors are removed, so a steady-state simulation is obtained; for the powered case, the propulsor is rotating providing an unsteady simu- lation. Even though the nonlinear free surface implemen- tation assumes a steady-state flow, it is assumed that the unsteadiness produced by the propulsor can be adequately In the following section, the free surface algo- nthm is presented, which includes the method for solving the kinematic free surface equation, the imposition of the hydrostatic pressure on the Navier-Stokes equations and the grid movement algorithm. Then, the results for the DTMB Model 5415 series simulations are presented. NAVIER-STOKES SOLUTION ALGORITHM Ihree different sets of equations are solved in the process of simulating the flow around surface vessels and sub- merged vessels near the water surface. The unsteady in- compressible Reynolds-averaged Navier-Stokes equations, which are presented here in Cartesian coordinates and in conservative form, are solved to determine the velocity and pressure within the computational domain. Either the Spalart-Allmaras or the q—w turbulence model is used to simulate the turbulent viscosity primarily within the boundary layer, and the kinematic free surface equation is solved to advance the free surface in time. GOVERNING EQUATIONS Assuming that gravity acts in the y-direction (i.e., the ver- tical direction), the incompressible Navier-Stokes equa- tions can be expressed in dimensional form, denoted with superscript *, as 0~* ~v* dw* dx* by* ~z* ~t* + `'x* (u*2 + p ) + ,'0* (u*v*) + ~ (u*w*~=~*V2u* dv* + ~ ~u*v*y + 00* (V*2 + p + 9Y +~*(v*w*~=~*V2v* oot* + oo* (u*w*) + o0* (v*w*) + Ad (W*2 + P ) = p*V2w* oz* ~ POO, (2) The variables in the preceding equation are normalized with respect to a characteristic length scale (L) and free stream values of velocity (UOO), density (pOO), and vis- cosity (,uOO). Thus, the Reynolds number is defined as Re = UooLI z/OO. Pressure is normalized with P = (P* + Poo9Y*—poO)/pOOU2, where P* is the local dimensional static pressure. Following (Chorin, 1967), an artificial time derivative term (63pa/5t, where Pa = P//~3) has been added to the continuity equation to cast the complete set of governing equations into a time-marching form. The

nondimensionalized equations can be written in integral form as The solution algorithm consists of the following basic steps: reconstruction of the solution states at the `, _ 1 control volume faces, evaluation of the flux integrals for /0 Qd;V+ ~ F ndA = /^ G. n~dA (3) each control volume, and the evolution of the solution in (fit JO JAQ Re JOQ each control volume in time. where n is the outward pointing unit normal to the con- trol volume V. The vector of dependent variables and the components of the inviscid and viscous flux vectors are given as Q=~v I IW ~ 3(~-at) F n = up + nap wE) + n:P o G · n = need + ray + nz~Z n~z + nary + nz~z news + news + nzrzz where ~ is the artificial compressibility parameter (typi- cally 15 in this work), it, v, and w are the Cartesian ve- locity components in the x, y, and z directions, and no, no, and n: are the components of the normalized control volume face vector. (3 is the velocity normal to a control volume face: Ed = not + nyv + nzw + al (7) where the grid speed al = —(Vent + Vyny + V:n:). Note that Vs = Vat + Van + V:k is the control volume face velocity. The viscous stresses given in Equation 6 are Tij = (A + fit) (0Xj 0~z ) where ,u and ,~b~ are the molecular and eddy viscosities, respectively. NUMERICAL APPROACH FOR NAVIER-STOKES EQUATIONS The baseline flow solver is a node-centered, finite volume, implicit scheme applied to general unstructured grids with nonsimplical elements. The flow variables are stored at the vertices and surface integrals are evaluated on the me- dian dual surrounding each of these vertices. The nonover- lapping control volumes formed by the median dual com- pletely cover the domain, and form a mesh that is dual to the elemental grid. Thus, a one-to-one mapping exists be- tween the edges of the original grid and the faces of the control volumes. RECONSTRUCTION A higher order spatial method is constructed by extrapo- `4~ rating the solution at the vertices to the faces of the sur- rounding control volume. The unweighted least squares method (solved via QR factorization (Anderson, 1994~) is used to compute the gradients at the vertices for the ex- trapolation. With these gradients known, the variables at (5) the interface are computed with a first order Taylor series as: Qf = Qo + VQo r `9y where r is a vector that extends from node 0 to the mid- `6' point of the edge associated with the control volume face In question. RESIDUAL EVALUATION The governing equations are discretized using a finite vol- ume technique; thus, the surface integrals in Equation 3 are approximated by a quadrature over the surface of the control volume of interest. So, the numerical discretiza- tion of the spatial terms associated with the control vol- ume surrounding vertex 0 results in Sqo + ~0 = 0 (10) where the spatial residual ~ contains all contributions from the discrete approximation to the inviscid and viscous terms (~ = Kiev + Evict. Also, the quantity q is defined as q = By Qd;V. SPATIAL RESIDUAL The evaluation of the discrete residual is performed sepa- rately for the inviscid and viscous terms given in Equation 3. The Roe scheme (Roe, 1981) is used to evaluate the in- viscid fluxes at each face of the control volume. The alge- braic flux vector is replaced by a numerical flux function, which depends on the reconstructed data on each side of the control volume face: ~ = 2 (IF f QI ~ + FTQR))—2 A (QR—QI ~ (1 1) ~ ~ ~ ~ where A = CAR-. The matrix R is a matrix constructed from the right eigenvectors of the flux Jacobian, and A is a diagonal matrix whose entries contain the absolute val- ues of the eigenvalues of the flux Jacobian. The eigensys- tem used in this work is based on that reported in (Tay- lor, 1991~. Note that A is evaluated with Roe-averaged

variables, which is simply the arithmetic average between left and right solution states in the case of incompressible flows. For general element grids, it is expedient to use only edge-local information to compute the viscous fluxes. This allows the evaluation of viscous fluxes on each face of the control volume without regard to the varying ele- ment types of the mesh. An algorithm in which no ele- ment information is used outside of metric computations is termed a "grid transparent" algorithm (Haselbacher, 1999~. To this end, the viscous fluxes are evaluated directly at each edge midpoint using separate approximations for the normal and tangential components of the gradient vector to construct the velocity derivatives (Marcum, 1997~. Us- ing a directional derivative along the edge to approximate the normal component of the gradient and the average of the nodal gradients to approximate the tangential com- ponent of the gradient leads to the following expression (Hyams, 2000b): VQij ~ VQ + (Qj - Qi - VQ As] -2 (12) The weighted least squares method is used to evaluate the nodal gradients in the preceding formula. TEMPORAL RESIDUAL After the spatial terms have been suitably discretized, the time derivative term appearing in Equation 10 must be ap- proximated. A general difference expression is available for this purpose (Beam, 1978) (Taylor, 1993), and is given as follows: /`qn = B1/Yt 0(Aqn) 1+62 fit (13) ~t ~ (qn) + 82 Con—1 where /`qn = qn+l _ qn. A first order accurate in time Euler implicit scheme is given by the choices b1 = 1, 62 = 0. Correspondingly, a second order time accurate Euler implicit scheme is given by B1 = 1, 82 = 1/2. Since 81 = 1 for both time discretizations used in this work, Equation 13 can be further simplified: In = /\t 63 (qn+l) + 82 Akin—1 (14) Using Equation 10 to replace the time derivative, Ax — 1+92 my 1 ~n+1 <15y At 1 + 82 where By the definition of q, one can write q = QV, where Q is ~0 + (Q + ) the volume averaged solution variable vector 1/V iso QdV. Then, the following two identities can be formed: /`qn—1 = An—1~\Qn—1 + Qn~`vn—1 (17) Inserting the above two identities into Equation 15, one arrives at the following expression: Vn+1 /\ Qn _ 82 in—1 /\ Qn—1 /\t Q [ At ] (18) + 1 ~pn+1 = 0 1+82 Now, one must consider the Geometric Conservation Law (GCL). This statement relates the rate of change of a phys- ical volume to the motion of the volume faces: Alit = ~ V · V~dV = ~ Vs ·ridA (19) According to Thomas and Lombard (1978) and later Janus (1989), the solution of the volume conservation equation must be performed in exactly the same manner as the flow equations to ensure that GCL is satisfied. This proce- dure ensures that spurious source terms caused by volume changes are eliminated. Using the same time differencing expression (Equation 14) to approximate Equation 19, vn_ 82 In—1 1 1+82 = 1 + ~ GIL (20) where GIL = I(0) VnO+i1 · ~oi+l' and \(0) is a list of nodes surrounding node 0. Note that the left hand side of the preceding equation is exactly the bracketed term in Equation 18. Replacing the bracketed term and rearrang- ing slightly gives the final form of the discretization of the time derivative: (1 + 6~2 )vn+1 is Qn—82Vn—1 ~ Qn—1 2\t QngRn+1 + jRn+1 = 0 (21) For incompressible flows, a divergence-free velocity field at the end of each Newton iteration is desirable. To this end, the contribution of the time derivative and GCL terms to the residual are removed for the continuity equation. TIME EVOLUTION A Newton iterative time evolution scheme is applied, which requires the solution of a sparse linear system at each non- linear subiteration: duct n+1 m _~On+l~m = U0 ~~Qn+l,m (22) ( 1 + 62 ) VOn+l i\ QOn—§2 VOn—1 i\ QOn—1 ~23' At /`qn = Vn+l /`Qn + Qn~vn (16) + Q ~O,GCL + ~0

Qn+1 m Qn+l,m+1 _ Qn+l,m. Now, ex- panding the terms and performing the required differen- tiations of ~ results in the following expression for New- ton's method: ~ (1 + 82 judo + (Qon+ltm _ QOn)—82vn—1~Qn—1 _ At + QO.~O,GCL + ~ Hoi ick(O) (1 + 82)Vo +1I L At + ~ `~j~t am . rlOni+1 1 ~ +1 - The initial guess to beg~n the iteration is /\Qn+l,° = 0. In the above formulas, Hoi represents the combined invis- cid and viscous flux vector from node O to node i. \~0) represents the set of neighbors for node 0, N~.(O) is the list of neighbors such that the node label earl < Eros, and ~u(O) is the list of neighbors such that the node label earl ~ egos BOUNDARY CONDITIONS . rl,r~+1 I = Viscous conditions are enforced by modifying the linear system such that no change is allowed in the velocity, and the pressure is driven according to the imbalance in the continuity equation in the boundary control volume (Anderson, 1994~. Farfield conditions are handled via a characteristic variable reconstruction; all boundary condi- tions are handled in an implicit fashion. The free surface boundary conditions will be presented below - in essence, ~~Hn+l,m -n+1 ~ at steady-state, the velocity is tangent to the free surface iEN~(O) ~ Oi no, ~Qu+ltm] (24' diver, ~ nfs = 0) Ad Ale pressure is P = An, where Y where /`QOn-1 = QOn _ QOn-l. For notational conve- nience, both the inviscid and viscous terms are collapsed into a single flux function H. Note that the iteration can be started by using an initial guess of Qn+1,0 = Qn. Also, performing only one iteration of Newton's method per time step (with 1st order time discretization and no GCL terms) is equivalent to a time linearization of the spatial terms only. However, writing the method in this framework is more general than a straightforward time linearization of the nonlinear terms. To solve the resulting linear system, a bidirec- tional Gauss-Seidel solution algorithm is adopted. Split- ting the matrix into diagonal, upper triangular, and lower triangular parts as [~] = [L] + [~] + tU], the Gauss-Seidel sweeps may be written as the following two-step process per subiteration (k is the linear subiterative index): tr + v~ ~Qk+2 + tu~ ~Qk = ~n+ltm (25) tv + Up ~Qk+1 + tr~ Auk+- = ~n+~,m (26) where it is understood that /`Q is evaluated at the rt + 1, m + 1 time level and Newton subiteration. The diago- nal, lower, and upper operators are defined as D ~1 + 82)VoI + L At dHo7li+l~m . rl,On+l ~ ic~O) Who ~ () `,Hn+l,m -n+i = ~ Oi · Hi (.) L= ~ iced (o) is the elevation of the free surface above the undisturbed waterline and the Froude number is Fr = Add. A symmetry plane boundary condition imposes V; · b = 0 for any arbitrary variable ¢. In addition, no flow is allowed through a symmetry surface; so, like a solid wall condition, ~ = 0. To simulate symmetry condi- tions, a layer of phantom entities that is a mirror image of the entities inside and connected to the symmetry plane is created; one element layer is mirrored. Thus, control vol- umes on the symmetry plane are closed and behave just as interior control volumes. Phantom nodes created by the mirroring process are updated in the same hierarchy as described for phantom nodes generated by the paralleliza- tion. Note that care must be taken to copy scalars, mirror vectors, and mirror tensors appropriately to the symmetry phantom nodes such that no fluxes are allowed through the plane. TURBULENCE MODEL Both the one-equation turbulence model of Spalart and Allmaras (1992~(Sheng, 1999) and the q—w turbulence model (Coakley, 1985) are available within the solver. Con- stants are taken from the version II q—w model given in (Coakley, 1985~. The interested reader is referred to these references. The diffusive terms in both turbulence models are discretized in the same manner as the viscous terms for the mean flow, and the convective terms are computed <27y via pure unwinding. Appropriate consideration is given to maintain positive operators in the formation of the Jaco- (28) bian matrix for the implicit solution of the transport equa- tionfs). The respective turbulence models are incoroo- . I . ] . 1 ~ ~ n+l~m nr,.+1 rated with the mean flow solution in a "loosely-coupled" Oi Em oi `.y (29) procedure;thatis,the core governing equations are solved

I ~ ~ ~ ~ ~ ~ in, rQ ~q At I ;::: lj';^IA ~ ~ ~ :: ~ ~ L ~ ~ J s~ ma : ~ : :: : ::~:__ r ~ t ~ em l ~ jut T;m" atnr. ~ ~ i:: :.: . :: : : : ~ : :~: ~ l ~ ~ ~~ v—_ ' ~ ~ ~ ~ Ha [Q. VQ, <* I Cycle over local subdomains~:~:: i: ~ I r- * ~ Lit, lit] ~~: : A? ~ T:~ ~ ?~: : ~ _ NewtonIteration : ~ :~:~ : ~~ Compute VQ em VQt Compute ~ Compute [A] ~ I Linear ~ subiterations: ~: : ~ ~ ~ :1 _- Matrix-Vector Muliiply - _ A <)t ~ ~ ~ :~: ~: : : ~ ~ ::~: ~ :: ~ : ~ ~ ::: :: :~: ~ 1 ~~SolYe~I\frbulence~:~Model : ::~: ,,,,B,,~ : :: : ~ [u pt]t ~ ~~ ~ ~~ ~ ~ ~ 1 ~ ~~ol~re:~::~nemat~c free Menace - ~ t _ _—[Yf~~Pfs] * sequential updating hierarchy t concurrent updating hierarchy Figure 1: Iteration hierarchy used for the parallel unstruc- tured solution algorithm first, then the turbulence model is solved independently. This procedure allows for easy interchange of the turbu- lence models. PARALLELIZATION For quick turnaround time in a design environment, it is essential to parallelize the flow solution algorithm. The present parallel unstructured viscous flow solver is based on a coarse-grained domain decomposition for concurrent solution within subdomains assigned to multiple proces- sors. The solution algorithm employs iterative solution of the implicit approximation, with the concurrent iteration hierarchy as shown in Figure 1. Also, domain decompo- sition takes place with each node in the domain uniquely mapped to a given task. The code employs MPI message passing for interprocessor communication. In general, the parallelization of an existing val- idated flow solver should satisfy several constraints. First and most important, the accuracy of the overall numerical scheme must not be compromised; i.e., the solution com- puted in parallel must have a one-to-one correspondence with the solution computed in serial mode. For the cur- rent numerical algorithm, this ability has been shown in (Hyams, 2000~. Also, the code must be efficient in its use of computational resources. This characteristic is mea- sured in terms of memory usage and scalability, as well as the fact that the parallel code should degenerate to the serial version if only one processor is available. Finally, the consequences of the inevitable domain decomposition should not seriously compromise the convergence rate of the iterative algorithm. NONLINEAR FREE SURFACE ALGORITHM A nonlinear free surface is obtained for a steady-state sim- ulation by solving the kinematic free surface boundary condition at each time level and imposing the hydrostatic pressure distribution based on the new free surface eleva- tion onto the free surface boundary within the mean flow. After several time steps, typically on the order of 200, the grid is moved to match the free surface elevations while conforming to the surfaces that intersect the free surface, such as the hull of a ship or the sail of a submarine; and as the loosely coupled interaction between the free surface, the Navier-Stokes equations and the grid movement al- gorithm converges, the solution approaches the nonlinear free surface solution. GOVERNING EQUATIONS In deriving the kinematic free surface boundary condition, the free surface ~(x, y, z, t) = 0 is considered as a mate- rial surface (i.e., no flow travels through the surface), so D71 = 0 (30) or lilt + (a- Vx) ,~,71 + (v- Vy) '971 + (w-Vz) ,0~71 = 0 (31) By making the assumption that the free surface 77 can be expressed as 71(X, y, z, t) = y—Y(x, z, t) = 0, this equa- tion reduces to lit + (a - V=) `~ - (v - Vy) + (w - Vz) ,, = 0 (32) At steady-state, the kinematic boundary condition becomes (u+at) (~,-1' Liz) 0 (33) (U+at)ootrifs =0 where or = )/1 + ( OoY ) + ~ BY ~ , so the flow is tangent to the free surface at steady-state. The pressure in the numerical simulation at the free surface is set based on the free surface elevation. For the original, dimensional equations, the pressure at the free surface is the atmospheric pressure POO- Using the nondimensionalization, p = P +Pp°°9u2-P°° with P* = POO and ~ = Y. the pressure is clearly set to P = An, where the Froude number is up . NUMERICAL APPROACH FOR FREE SURFACE The kinematic equation for the linear free surface is solved via a Galerkin approach. In this approach, an algebraic equation is obtained for each node j by multiplying the

governing equation by a weight function ~j(x,z), ap- proximating the terms in the governing equation by using linear and bilinear interpolating functions Minx, z) and in- tegrating over the computational domain, or ~Q~i(X~Z) (at + (U-V~) o~ - (V_Vy) ~ \ +(W-Vz),g )=0 (34) where the interpolated free surface elevation is Y(x, z) = Pi Yi (i (x, z) and the velocities are of the form USA, z) = Ui~i(X~Z). For the triangular elements, the resulting equa- tions are integrated either exactly or using 7 point Gauss quadrature, and for the quadrilateral elements, the inte- grals are solved either using 3 or 6 point Gauss quadra- ture in both directions or by decomposing the quadrilat- erals into 4 overlapping triangles and evaluating them as triangles. The examples given herein use 7 point Gauss quadrature for the triangular elements and 6 point Gauss quadrature for the quadrilateral elements. A first-order backward time discretization of Equation 34 is used for the temporal derivative, where CITY = Y lit Y and the spatial terms evaluated at time level r' + 1. This dis- cretization results in a linear algebraic system of the form ~fs(Yn+t, Y~) = 0. Since these equations are linear in the unknown variable Yin+i at each time level (i.e., the ve- locities are frozen), the Jacobian matrix is calculated only once per time level. A Gauss-Seidel iterative solver simi- lar to the one used to solve the discretized Navier-Stokes equations is used to identify the solution of the kinematic free surface equation by solving of the viscous surface is on the order 10-6 for model scale simulations and 10-8 for full scale simulations. FREE SURFACE BOUNDARY CONDITION At the free surface boundary, the pressure is set to the hy- drostatic pressure of P = ma. Determining the velocity to impose on the boundary is more difficult. If the ve- locity were forced to be tangent to the free surface (i.e., (n + at) · Ifs = 0), then the free surface equation would reduce to MY = 0, and the free surface could not evolve. Thus, this tangency condition is relaxed to allow the free surface to evolve. This boundary condition is implemented using a characteristic variable boundary condition, similar to the derivation used for the farfield boundary condition. For hyperbolic systems, flow information travels along characteristics; and for three-dimensional incompressible flow, three characteristics originate upstream, and one orig- inates downstream. If the flow is traveling out of the do- main (i.e. (u + at) · n > 0), then only one characteristic needs to be specified from the outside, which is derived from the hydrostatic pressure. If the flow is traveling into the domain, three characteristics need to be specified from the outside. The characteristic variables are derived from considering the inviscid fluxes across the the boundary face, via ~Q + VF rl = 0 (36) where r' is the normal to the grid for each boundary face. The inviscid flux term can be rewritten as VF = L3Q VQ = RAR-I VQ (37) where A is a diagonal matrix with the eigenvalues 0, 0, `~ (3 + c and (3—c, where (3 = uric + vr~y + wry + at and n+t,m _ ~ ~yn+i,m yin ~35> ~ 2 ~yn+i Y —— fey ~ J ~ J c=~(~)_ I) +~. PremultiplyingbyR-t, where yn+i,m+i = yn+i,m + /`yn+i,m and yn+i = -~0Q AR-IV —0 38 yn+itm+i when ~~Ayn+~,m~ < tolerance. R ~ + Q ~— Within the viscous boundary layer, special care is needed in solving the kinematic free surface equation. On viscous surfaces, the flow velocity plus the grid ve- locity (i.e., u + all is set to zero, which prevents the free surface from moving at the viscous surface. Physically, however, surface tension effects force the free surface to rise and fall along the viscous surfaces. Since the solver does not simulate surface tension, another method must be employed to move the free surface near the viscous sur- faces. Following the method used in the structured solver, the free surface at a certain distance from the viscous sur- face is extended at a constant height to the viscous surface. This distance is typically on the order of 2 x 10-4 times the characteristic length L, whereas the point spacing off Freezing the matrix Ro ~ = R-t f~Qin`), ARE Q +~V (Ro 1Q) ri = 0 (39' Using the definition W(`Q`) = Ro-iQ' Equation 36 decou- ples into four equations for the four characteristic vari- ables W as 0~9~ +AVW ri = 0 (40) Three of the characteristics are obtained from upstream values, and one of the characteristics is obtained from downstream values. In the case of a boundary face, there are a set of characteristics associated with the inside flow

I. Wln W = win w4ut characteristic variables are chosen as we) and a set of characteristics associated with the into the volume grid, the use of a linear spring analogy outside flow W°Ut (Q°Ut). The flow variables on the bound- where each edge is a spring whose stiffness is determined ary face are determined from the appropriate characteris- by the length of the edge, a torsional/linear spring anal- tic variables via Q = RoW. For flow traveling out of the ogy where the angle between the edges affects the stiff- domain, the characteristic variables are chosen to be ness of the springs in the mesh, or solving the linear elas- i~ ticity equations to propagate the perturbations within the wag grid. The spring analogy is computationally efficient but W = w2n <41y is not robust. Both the spring analogy and the Laplacian wout solver generate elements with negative volumes for mod- 4 crate amounts of movement, on the order of the size of 3 because the fourth eigenvalue (3—c is negative, while to 4 elements. Solving the linear elasticity equations to the other eigenvalues are positive. For flow traveling into move the grid is robust, but the computational cost asso- the domain, only the third eigenvalue is positive, so the elated with this method is the primary drawback. Thus, following the work of (Farhat, 1998) and extending it to three dimensional grids, the torsional/linear spring analogy has been developed and applied to the prob- `42' lem of moving the grid to match the linear free surface el- evations. For tetrahedral grids, this method is quite robust, allowing severe distortions on the surface while maintain- ing positive volume elements. In practice this algorithm For farfield boundary conditions the characteristic vari- . .' ' allows up to approximately 80% compression of elements. ables associated with the outside flow are taken as the . . Within the boundary layer, (led, near the viscous surfaces), freestream variables. For free surface boundary cond~- .. . . the grid moves the same as the points on the boundary. tons the characteristic variables for the outside flow are ~ This grid movement method Is computat~onally costly but determined from the hydrostatic pressure and the only avail- prov~des excellent robustness for the free surface solver. able velocity information which is the velocity on the ~n- ~ The linear so arms method for moving an unstruc- side. However this velocity is modified by removing a . . . ~ ~ ' tured arid Is presented In (Bating, 1989) where each edge portion of the velocity component that is not tangent to . ~ . . . . . . In the cad Is replaced by a linear spring whose stiffness Is the grid, via I-WO?lt wout wont ~ = (A + at)—by ((n + at) n) (n + at) (43) where :_ ~—(Nat) ~ if —1 < (~+at) ~ < 0 t1 if (Nat) ~ < -1 (44) When y is 1, the velocity imposed from the outside is tan- gent to the grid (i.e., clout ~ = 0~. When the velocity was not modified or constrained in some fashion, the free surface algorithm became unstable for flow into the do- main due to the inconsistency of using downstream infor- mation. This modification provides enough control over the velocity to maintain stability of the algorithm, and if the grid is allowed to move to match the free surface, then this modification ensures that the flow will be tangent to the free surface at convergence. GRID MOVEMENT ALGORITHM After several time steps, the grid is moved to match the free surface while conforming to any solid surfaces inter- secting the free surface, with displacements on the surface being propagated into the volume grid. Several methods are available for this grid movement, including the use of a Laplacian solver to propagate the surface perturbations inversely proportional to the length of the edge. Thus, for the edge connecting nodes i and j, the stiffness kid of the . . spring IS ( (Xi—X j ) 2 + (hi - y j ) 2 ) p/2 = ],[ p (45) where the coordinates of nodes i and j are (xi, pi) and (Xj ~ yj) and p is a predetermined coefficient, usually 1 or 2, and lid = \/((xi—Xj)2 + (Pi—yj)2). Given a set of nodal displacements or forces acting on the boundary of the computational grid, the following equations for the interior displacements are solved iteratively until all the ~ . . . . forces are In equip. ~ Pram i ~ k i ti ~yn+1 = E; kid EVE i 3} 1 (46) where i is summed over all edges connected to node j. This method is easily extended to three dimensions and is computationally efficient requiring only a few Jacobi iterations to achieve an acceptable level of accuracy. Sev- eral researchers have used this method to move nodes in unstructured grids for simulating flows around objects in relative motion (Singh, 1995) and for design optimization (Anderson, 1997).

./~a (xi,yi) Xj-Xi ~ (Xj ,yj ) Yj-Yi 2 Figure 2: A Displacement in y-direction results in dis- placements in x and y directions Technically speaking, these equations do not sim- ulate the behavior of a network of springs because there is no interaction between the x and y coordinates. A dis- placement in one coordinate will not influence the loca- tion in the other coordinate, as would be the case for a network of springs. An example of this coupled inter- play is shown in Figure 2, where a displacement in the y-direction is applied to node j which results in a displace- ment in both directions for node i. To simulate the behav- ior of a network of springs, the following set of equations for each spring must be summed over all springs: forces in the x-direction: kij [~Xi—Axj) cos2 ~ + (~\yi—/\yj) COS Ol sin Ol] <47> forces in the y-direction: kij [taxi—A~j) cos ~ sin + (phi—Ayj) sin2 a] For either representation, however, this method often fails for complicated geometries and for large changes in the boundary, because negative areas are generated when a node crosses over an edge in the grid. The creation of negative areas is illustrated in Figure 3, where node 1 is pushed downwards and node 4 is forced to cross the edge between nodes 2 and 3. The reason for this failure is that the stiffness in the linear spring method prevents two nodes from colliding but does not prevent a node from crossing over an edge. As two nodes get closer to- gether, the stiffness increases without bound preventing the collision; but there is no mechanism to prevent a node from crossing an edge, because these crossovers can occur without the stiffness increasing without bound. To provide a more robust movement algorithm for unstructured grid, Farhat (1998) developed an algo- rithm to prevent a node from crossing an edge by using torsional springs around each node, as shown in Figure 4. ~ ~L~ ~\~ 3 2 3 2 4 3 1 Figure 3: Negative Areas Produced by Linear Spring Method k Linear Spring :~ ~ \ ~ \ ~)X/\e ilk \ ire ~ . ~_J Torsional Spring Linear Spring Figure 4: Placement of Linear and Torsion Springs The stiffness of the torsion spring Cijk is inversely pro- portional to the sine of the angle, so that the stiffness grows without bound as the angle decreases towards zero or increases towards 180°, or Cijk = j 2 Jo,, (48) where 83jk iS the angle centered at node i formed from the edges ij and ik. Thus, as a node moves towards an edge, the angle goes towards zero, and the stiffness of the torsion spring grows. The sine of the angle is squared to prevent a negative stiffness. In (Farhat, 1998), the equations for the torsional spring methods are derived for two-dimensional grids, and an iterative solution method is presented. In the present work, these methods have been extended into three di- mensions, which provides an adequate level of robustness for the grid movement algorithm. The derivations of the three-dimensional torsional spring equations are beyond the scope of this paper and will be presented elsewhere. EXAMPLES The two examples given in this paper deal with the fully- appended DTMB Model 5415 series hull, one without propulsors and the other with propulsors rotating. Ex- perimental data for comparison was provided by Ratcliffe

at the Naval Surface Warfare Center Carderock Division (NSWCCD) (Ratcliffe, 1990, 2000), which consisted pri- marily of the profile of the free surface intersection with the hull and free surface elevations in the stern of the hull. The Model 5415 series hull has a transom stern, which produces a sizeable "rooster tail" in the water surface be- hind the stern. The parameters for the various runs are given in Table 1. r Parameter | Value l ~ I 1 1 Speed TOO (m/sec) 2.06 Ship Length L (m) 5.72 Propulsor RPM 436 Advance Ratio (Rev. per L) 20.177 Froude Number (Fr) 0.28 Reynolds Number (Re) 11.655 Million Omega 126.59 red/second At for 1.5 degrees 2.068089 x 10-4 Table 1. Flow parameters for Model 5415 Series simulations. The grids used within the following simulations were gen- erated by AFLR3 (Marcum, 1998) and consist of prisms in the boundary layer generated from the triangles on the viscous surfaces, tetrahedra and pyramids in the transition region at the edge of the boundary layer grid and tetrahe- dra outside of the boundary layer. On the free surface grid, quadrilaterals exist within the boundary layer region and triangles exist outside of the boundary layer. The use of tightly packed prisms within the boundary layer allows for good resolution of the boundary layer effects and typ- ically generates y+ values on the order of 1.0. Information about the size of the grids used in these simulations is provided in Table 2. For the powered and unpowered simulations, the same grid is used except for the regions around the two propulsors. For the un- powered grid, the propulsors were replaced by hubs and were not rotated. For these simulations, the grid included both sides of the hulls, so that maneuvers using these grids could be studied in the future. The hull, struts, rudders and propulsors were symmetric about the centerline, but because of the use of unstructured grids, the grid was not symmetric about the centerline. These simulations were run locally on a LINUX-based cluster running MPI-Pro. I l Unpowered l Powered l Hull Profiles for Model 5415 Hull 0.018 0.015 0.012 0.009 0.006 0.003 o -0.003 -0.006 0 0.2 Fr=0.28, Re=11.355 M ' I ' I ' 1 ' 1 ' Barehull (Structured) - Unpowered - Powered 0 DTMB Experimental Data 0.4 0.6 0.8 1 Distance from Bow of Ship Figure 5: Free Surface / Hull Intersection which the simulation was run at a constant time step, so that a time accurate simulation could be achieved for the powered case. To ensure a valid comparison, both simu- lations were run at a constant time step. For the powered case, the images and plots use instantaneous data rather than time-averaged data. UNPOWERED FULLY-APPENDED 5415 The unpowered fully-appended 5415 simulation used the same grid as the powered fully-appended 5415 simulation with the exception that the grid around the propulsor/hub was changed and did not rotate. Hence, the free surface grid was identical and the forces acting on the free surface near the bow should be almost identical between the sim- ulations. The profile of the intersection between the free surface and hull is given in Figure 5. This profile com- pares the results for the powered and unpowered simula- tions with the experimental results from the Naval Surface Warfare Center Carderock Division runs (Ratcliffe, 1990, 2000) and with the computational results of the structured flow solver UNCLE (Beddhu, 1999, 2000~. The profile for the powered and unpowered cases agree with each other nicely until approximately y/L = 0.6 at which point the unsteady influence of the propulsors are felt. The bow wave in the unstructured grids are higher than for the structured grid and seem to agree more with the bow wave elevation of the experiment. The drop behind the bow wave is more severe for the unstructured simulations Nodes 4251444 5437074 - this difference may be caused by a lack of refinement in Tetrahedra 8195672 9912263 this region. Pyramids 23592 35101 Figure 6 shows the bow wave for the unpowered Prisms 5446868 7160111 simulation head on, as well as the bulbous bow underneath Processors 32 42 the free surface. In Figure 7, the surface grid near the bow Table 2: It formation abo at Grids. is shown. This region undergoes the largest displacement due to the changes in the free surface, and yet the grid In each case, the simulation was run with a constant CFL movement algorithm is able to produce nice quality in the until the gross features in the flow field developed, after

Figure 6: Bow Wave for Model 5415 Figure 7: Grid Moved to Match Free Surface Elevation near Bow distorted grid. The hullform was symmetric about the center- line, but due to the use of unstructured grids, neither the volume nor the surface grids were symmetric. This asym- metry in the grid produces an asymmetry in the discretiza- tion error within the code, but the free surface, even in the stern region, regains symmetry as the solution approaches steady-state. In Figures 8 and 9, contour plots of the free surface elevations around the unpowered 5415 are pre- sented showing the Kelvin wake pattern originating from the bow and the "rooster tail" in the stern region. For the unpowered case, the solution is nearly symmetric. Due to the coarseness of the grid away from the body, the free surface dissipates approximately one body length behind the hull. By extending the resolution on the free surface, the free surface waves could be carried further. The free ..f ~ , Figure 8: Contours of Free Surface for Unpowered 5415 Figure 9: Contours of Free Surface for Stern Region of Unpowered 5415

Unpowered Model 5415 Stern Topology O.05 -D.05 —1.05 1.1 X,L y 0.~e 0.~07 0.006 0.005 O.~04 0.003 0.002 0.00 1 .~D 1 .002 .003 Figure 10: Comparison with Experimental Data for Un- powered 5415 surface in the stern region was also compared Wltn tne ex- perimental data from Ratcliffe at NSWCCD and is shown in 10. The agreement for the unpowered case is remark- able, with similar features and elevations. There is a no- ticeable shift forward for the computational results. Finally, the free surface-stern intersection for the powered and unpowered simulations are given in Figure 11, as well as the profile for the structured barehull simu- lation. The structured barehull simulation was run using a symmetry plane, so its profile was mirrored about z = 0. The unpowered simulation is almost symmetric about the stern, whereas the unsteadiness in the powered simulation has created some asymmetry (this asymmetry is probably a result of not being as highly converged as the unpowered case). For the powered simulation, this stern profile is an instantaneous cut, and is not as converged as the unpow- ered simulation. The profiles for the simulations differ greatly. The reason for this difference lies primarily with the difference in the flow field behind the rudders. For the barehull, the struts, shafts and rudder have been removed, so the how has not been slowed nor straightened, thus the flow passing the stern is faster for the barehull case than for the unpowered case. For the unpowered case, the flow has been slowed by these appendages and straight- ened by the rudder, so that the flow is slower, and the free surface is not sucked down. And for the powered case, the propulsors have accelerated the flow past the rudders, which sucks the free surface down. In Figure 12, the ver- tical velocity on the free surface and for a cutting plane just aft of the stern are presented from a view beneath the hull, showing that the velocity is relatively uniform be- tween the rudders and the vertical velocity is nearly uni- . .. . . . . · .. Free Surface Stern Intersection 0.002 ,, / ,, -0.002 I A v -0.004 -0.04 -0.02 i.\ ,j ~ ~ , ~ —,j ,~ ~ if ,' -it—·_,,, ., " ll Powered Unpowered ---- Barehull I' . 0 0.02 z/L Figure 11: Free Surface / Stern Intersection 0.04 Figure 12: Vertical Velocity for the Unpowered 5415. form along the free surface behind the stern. In Figure 13, the instantaneous vertical velocity for the powered case is shown from the same location and on the same scale. The swirl in the velocity field generated by the propulsors are clearly seen and indicates that the propulsors are rotating outboard over the top. The velocities are much more ex- treme at the free surface behind the stern, which results in the drop in the free surface elevation along the stern. POWERED FULLY-APPENDED 5415 The powered fully-appended 5415 simulation was per- formed using the same parameters as the unpowered case, with the exception that the powered case rotated the ac- tual propulsors used in the experiment. As a result, a pe- riodic solution is expected for the powered case. Figure 14 shows the drag on the port propulsor as well as the to- tal drag on the hull, struts, shafts and rudders and on the

0.0008 ~ 0.0004 c, n I I l ~ / / / Figure 13: Vertical Velocity for the Powered 5415 (Instan- Figure 15: Contours of Free Surface for Stern Region of taneous). Powered 5415 0.0012 — Drag (Thrust) for Powered 5415 Nondimensional - Powered Fully-Appended 5415 Stern Topology 51 52 53 54 Rotation Count Figure 14: Drag (Thrust) on Powered 5415. Overall 5415. The drag on the propulsors is negative be- cause the propulsors are introducing thrust into the simu- lation. The drag is clearly periodic, approximately repeat- ing every 1/5 of the propulsor rotation, since the propul- sors had 5 blades. The overall drag is reducing slightly in- dicating that a highly converged solution has not yet been achieved. A contour plot of the free surface in the stern re- gion for the powered case is shown in Figure 15 and shows that the flow field is nearly symmetric. Finally, the stern region for the powered case is compared against the exper- imental results from NSWCCD in Figure 16. The exper- imental results are time-averaged, whereas the computa- tional results are instantaneous. The height of the "rooster tail" for the computational results is higher than for the experimental results. This difference may be caused by V3 · 0.0080 0.0072 0.0064 0.0056 0.0049 0.004 1 0.0033 0.0025 0.0017 0.0009 0.000 1 -0.0006 -0.00 1 4 -0.0022 _ -0.0030 55 56 1.D5 1.1 x`,Ll .15 1.2 Figure 16: Comparison with Experimental Data for Pow- ered 5415 the assumption in the nonlinear free surface algorithm that the free surface is steady-state when the flow is clearly un- steady. CONCLUSIONS A robust, nonlinear free surface algorithm coupled with an incompressible Navier-Stokes solver on three-dimensional unstructured grids has been presented. The underlying un- structured solver has been developed to simulate turbulent boundary layer effects and is capable of performing rud- der and propulsor induced maneuvering. The addition of nonlinear free surface capabilities is a necessary step to- wards the simulation of maneuvering surface vessels. The nonlinear free surface algorithm is an ex- tension of the algorithm used in the structured code UN-

CLE developed by researchers at Mississippi State Uni- versity. The kinematic free surface equation is solved at each time level, and the hydrostatic pressure is imposed on the Navier-Stokes solver via a characteristic variable boundary condition. After several iterations, the grid is moved to match the free surface while conforming to the underlying geometry, using a three dimensional extension of the torsional spring method (Farhat, 1998~. By allow- ing the grid to move, a nonlinear free surface solution is obtained. In this paper, the free surface around the fully- appended DTMB Model 5415 series hull with and with- out rotating propulsors has been simulated. The hull pro- files for both simulations upstream of the propulsors is in close agreement with each other and adequate agreement with the experimental data and with computational results using the structured solver UNCLE. The contours of the free surface elevations for both simulations show that the free surface simulation is symmetric even though the grid was not symmetric. And the flow in the stern region for both geometries is symmetric and agrees well with the ex- perimental data generated by Ratcliffe at NSWCCD. Dif- ferences in the free surface at the stern can be attributed to differences in the velocity field between the propul- sors and between the rudders. More work is necessary in demonstrating that the powered simulation does have a periodic free surface, while the current work has shown that the drag on the geometry is nearly periodic. Future work includes performing these same sim- ulations at full scale Reynolds number, analyzing the ef- fects of stern flaps on the Model 5415 series hull forms and Series 60 CB = 0.6 ships, simulating flow around unconventional hull forms and performing prescribed and propulsor/rudder induced maneuvers. ACKNOWLEDGMENTS This research was sponsored by the Office of Naval Re- search with Dr. L. Patrick Purtell grant monitor. This sup- port is gratefully acknowledged. Special thanks to Toby Ratcliffe at Naval Surface Warfare Center Carderock Di- vision for providing the data files for the experimental data and to Stephen Nichols at the Computational Sim- ulation and Design Center for many hours of conversa- tions about free surface issues, especially pertaining to the structured code UNCLE. References [1] Anderson, W. K., and Bonhaus, D. L., "An implicit upwind algorithm for computing turbulent flows on unstructured grids," Computers in Fluids, Vol. 23, No. 1, Jan. 1994, pp. 1-21. [2] Anderson, W. K., and Venkatakrishnan, V., "Aero- dynamic Design Optimization on Unstructured Grids with a Continuous Adjoint Formulation," Proceedings of the 35th AIAA Aerospace Sciences Meeting end exhibit, AIAA Paper 97-0643, Jan. 1997. [3] Batina, J. T., "Unsteady Euler Airfoil So- lutions Using Unstructured Dynamic Meshes," 27th AIAA Aerospace Sciences Meeting AIAA Pa- per 89-0115, Jan. 1989. [4] Beam, R. M., and Warming, R. F., "An implicit factored scheme for the compressible Navier-Stokes equations," AIAA Journal, Vol. 16. No. 4, April 1978, pp. 393-402. [5] Beddhu, M., Taylor, L. K., and Whitfield, D. L., "A Time Accurate Calculation Pro- cedure for Flows with a Free Surface Using a Modified Artificial Compressibility Formula- tion," Applied Mathematics and Computation, Vol. 65, 1994, pp. 33-48. [6] Beddhu, M., Jiang, M. Y., Taylor, L. K., and Whit- field, D. L., "Computation of Steady and Un- steady Flows with a Free Surface Around the Wigley Hull," Applied Mathematics and Computation, Vol. 89, 1998,pp.67-84. [7] Beddhue, M., Jiang, M. Y., Whitfield, D. L., Taylor, L., K., and Arabshahi, A., "CFD Validation of the Free Surface Flow Around DTMB Model 5415 Us- ing Reynolds Averaged Navier-Stokes Equations," Proceedings of the Third Osaka Colloquium on Ad- vanced CFD applications to Ship Flow and Hull Form Design, 1998. [8] Beddhu, M., Nichols, S., Jiang, M. Y., Sheng, C., Whitfield, D. L., and Taylor, L. K., "Comparison of EFD and CFD Results of the Free Surface Flow Field about the Series 60 CB = 0.6 Ship," Proceedings of the Twenty-Fifth American Towing Tank Conference, 1998. [9] Beddhu, M., Jiang, J. Y., Whitfield, D. L., and Taylor, L. K., "Computation of the Wet- ted Transom Stern Flow over Model 5415," Proceedings of the Seventh International Conference on Numerical Ship Hydrodynamics 1999. [10] Beddhu, M., Pankajakshan, R., Jiang, M. Y., Taylor, L., Remotigue, M. G., Briley, W. R., and Whitfield, D. L., "Computation and Evaluation of CFD Results for Practical Ship Hull Forms," Proceedings of Gothenburg 2000: A Workshop on Numerical Ship Hydrodynamics, 2000.

[11] Burg, C. O. E., Sreenivas, K., Hyams, D., and [20] Mitchell, B., "Unstructured Nonlinear Free Surface Solutions: Validation and Verification" Proceedings of the 32nd AIAA Fluid Dynamic Conference and Exhibit AIAA Paper 02-2977 June , , 2002. [12] Chorin, A. J., "A numerical method for solving incompressible viscous flow problems," Journal of Computational Physics, Vol. 2, 1967, pp. 12-26. [13] Coakley, T. J., Hsieh, T., "A comparison be- tween implicit and hybrid methods for the calculation of steady and unsteady inlet flows." Proceedings of the 21st AIAA/SAE/ASME/ASEE Joint Propulsion Conference, AIAA Paper 85- 1125, July 1985. [14] Farhat, C., Degand, C., Koobus, B., and Le- sionne, M., "Torsional Springs for Two- Dimensional Unstructured Fluid Meshes," Computer methods in Applied Mechanics and Engin- [24] sing, Vol. 163, 1998, pp. 231 -245. [15] Haselbacher, A., McGuirk, J. J., and Page, G. J., "Fi- nite volume discretization aspects for viscous flows on mixed unstructured grids," AIAA Journal, Vol. [25] 37, No. 2, Feb. 1999, pp. 177-184. [16] Hyams, D. G., Sreenivas, K., Sheng, C., Bri- ley, W. R., Marcum, D. L., and Whitfield, D. L. "An investigation of parallel implicit solution algorithms for incompressible flows on multielement unstructured topologies," Proceedings of the 38th AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper 2000-0271, January 2000. [17] Hyams, D. G., "An Investigation of Parallel Implicit Solution Algorithms for Incompressible Flows on Unstructured Topologies," PhD thesis, Mississippi [27] State University, May 2000. [18] Hyams, D. G., Sreenivas, K., Sheng, C., Nichols, S., Taylor, L. K., Briely, W. R., and Whitfield, D. L., "An Unstructured Multielement Solution Algorithm for Com- plex Geometry Hydrodynamic Simulations," Proceedings of the Twenty-Third Symposium on Naval Hydrodynamics, September, 2000. [19] Janus, J. M. " Advanced 3-D CFD Algorithm for Turbomachinery," PhD thesis, Mississippi State University, May 1989. Marcum., D. L. Private communications, Engineer- ing Research Center for Computational Field Simu- lation, Mississippi State, MS, December 1997. [21] Marcum, D. L. "Unstructured Grid Generation Using Automatic Point Insertion and Local Reconnection," The Handbook of Grid Generation, edited by J.F. Thompson, B. Soni and N.P. Weatherill, CRC Press, Boca Raton, 1998, Section 18. t22] Ratcliffe, T. J., and Lindenmuth, W. T., "Kelvin Wake Measurements Obtained for Five Surface Ship Models," DTRC-89/038, January, 1990, David Tay- lor Model Basin. [23] Ratcliffe, T. J., "An Experimental and Compu- tational Study of the Effects of Propulsion on the Free-Surface Flow Astern of Model 5415," Proceedings of the Twenty-third Symposium on Naval Hydrodynamics, September 2000. Roe., P. L., "Approximate Riemann solvers, parameter vectors, and difference schemes," Journal of Computational Physics, Vol. 43, pp. 357-372, 1981. Sheng, C. Hyams, D., Sreenivas, K., Gaither, A., Marcu, D., Whitfield, D., and Anderson, W. K., "Three-dimensional incompressible Navier-Stokes how computations about complete configurations using a multiblock unstructured grid approach," Proceedings of the 37th AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper 99-0778, January 1999. [26] Singh, K. P., Newman, J. C., and Baysal, O., "Dy- namic Unstructured Method for Flows Past Multiple Objects in Relative Motion," AIAA Journal, Vol.33, No. 4, April 1995, pp. 641-649. Spalart, P. R., and Allmaras, S. R., "A one-equation turbulence model for aerodynamic flows," AIAA Paper 92-0439, 1992. [28] Taylor, L., K., "Unsteady Three-Dimensional In- compressible Algorithm based on Artificial Com- pressibility," PhD thesis, Mississippi State Univer- sity, 1991. [291 Taylor, L. K., Busby, J. A., Jiang, M. Y., Arab- shahi, A., Sreenivas, K., and Whitfield, D. L., "Time accurate incompressible Navier-Stokes simulation of the fiapping foil experiment." Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics, August 1993.

[30] Thomas, P. D., and Lombard, C. K., "Geometric conservation law and its application to flow compu- tations on moving grids" AIAA Journal, Vol. 17, No. 10, 1978, pp. 1030-1037.

DISCUSSION L. J. Doctors The University of New South Wales, Australia Thank you for a fascinating exposition of a very challenging and practical problem. I would like to ask you how you handle the low-speed flow calculations in the neighborhood of the transom stern. In this region, the water level on the surface of the transom must drop in a complicated fashion as the speed of the vessels increases, until some point is reached where the flow will separate clearly along the girth of the transom. AUTHORS' REPLY We choose a Froude number where the transom remained wetted. Thus, our initial "ridding efforts started with a planar free surface at the original undisturbed water line. Then, as the flow develops, the free surface sets up and the grid is moved to match the free surface. Around the corner of the transom, the free surface drops quite rapidly. If the grid were to match the free surface, the points on the transom would need to be compressed substantially. In this region, we relax the condition that the grid match the free surface, to enable us to have a robust grid movement/ nonlinear free surface algorithm. If we chose a Froude number which results in a dry transom, we would need to use a different initial grid. If we started with a planar free surface, as above, then we would need to compress the points along the transom down towards the bottom of the transom. Since we would not be able to compress them completely to an edge, the grid movement algorithm would fail and we would not be able to simulate this flow. However, by starting with a nonplanar free surface grid, we hope to be able to simulate the flow around a dry transom, as well as for other similar geometries. DISCUSSION Jacquin Erwan Bassin d'Essais des Carenes, France I congratulate the authors for the quality of their work, and for their presentation. Calculations using unstructured mesh, with nonlinear free surface conditions with rotating propulsors is a hard challenge, and the qualitative results you have shown seems very realistic. Nevertheless, no results are published on quantitative quantities that are available on this hull form in the literature (pressure and frictional resistance, boundary layer velocities, ...~. Did you make comparisions between your results and experimental data? Thank you for your answer. AUTHORS' REPLY You have clearly pointed out a weakness in our paper, in that we do not compare with experimental quantities such as drag or velocity profiles. The goal of the work that was presented was to develop efficient, robust and stable algorithms for the nonlinear free surface capabilities of our unstructured flow solver. In the future, we plan on analyzing the ability of our code to predict the forces on the hull form and the propulsors. Currently, our unstructured flow solver does not capture the hook feature in the wake region directly upstream of the propulsors. Our structured flow solver is able to capture this feature reasonably well. This failure is of great concern to us and we are seeking to remedy this. We also have concerns about the number of nodes required to achieve a grid independent solution. We have run the model 5415 series hull with several different grids and have noticed some significant differences in the flow field, along the hull behind the bow wave and in the stern region. Thus, we expect that we will need to run with significantly more points that we are currently using. Concerning the velocity distributions, we would not expect to match the free surface along the hull, if the velocity distribution near the hull was not accurate. And as we move away from the hull, which means that the grid is becoming less refined, we would expect that the velocity distribution would be smeared as is the free surface. By adding more points in these regions, we would probably maintain the velocity distribution better. Finally, the free surface in the transom region is highly dependent on Reynolds number which suggests that the free surface is dependent on the turbulent features in the flow. To adequately capture the flow features in this region, we will probably need to convect the turbulent features

upstream of the transom better, which includes simulating the hook feature in the nominal wake. As a predictive tool for forces and moments, we have much work ahead of us, before we and others will be confident of the results. DISCUSSION Toby J. Ratcliffe Naval Surface Warfare Center, Carderock, USA The authors have presented an impressive set of results from the Unstructured Grid code developed at MSU. Clearly, these results show much improvement in predicting the surface wave elevations behind the transom of Model 5415, in both the unpowered and powered conditions, over the structured grid solutions. At the Froude number of 0.28 for which these comparisons are made, the transom of the model is partially wetted. At higher speeds, the stern becomes completely dry. This dry transom condition has presented problems for modelling in a "traditional" structured grid code. Would the unstuctured RAN S code, detailed in this paper, with the non-linear free surface capabilty, be better suited toward solving the dry transom condition? AUTHORS' REPLY We have not run cases with a dry transom stern using either our structured or unstructured RANS solvers. The main issue for both codes is how to build the grid to satisfy the dry transom condition, and for the unstructured case, this issue is complicated with the issue of how to terminate the viscous packing off the viscous surfaces. For the wetted case, the viscous packing can terminate into the water surface without any special considerations, as shown in the diagram labeled Wet Transom Case. However, for the dry transom case, if the viscous packing were to terminate into the water surface, the grid quality at the transom corner would be quite poor. Other options for terminating the viscous packing have disadvantages as well. Hence, we would suggest wrapping the grid around the transom for one element, so as to allow the viscous packing to terminate cleanly into the water surface, as shown in the diagram labeled Dry Transom Case. By doing this, no special considerations for either the free surface solver or the grid movement algorithm would be required, and the points in the free surface along the transom stern would still be allowed to move up and down to a small extent and would then be a better indication of whether the numerics suggest that the transom is truly dry. Water Surface ,, ......... at Transom Hull Boy Wetted Transom Case r Transom | Hull Bo¢t~ DISCUSSION Dry Transom C:ase Michel Visonneau Ecole Centrale de Nantes, France / J water Surface Concerning the 3D regrinding tool, do you use the three dimensional generalization purposed by Farhat recently or did you build your own regridding tool? AUTHORS' REPLY In the fall of 1999, Dr. Farhat visited our university and shared with us his work in 2D torsional springs. We implemented his 2D algorithm and achieved similar excellent results. After several unsuccessful attempts to extend the method to 3D, we developed an algorithm which performs well in 3D. This algorithm was developed independently from Dr. Farhat during the winter of 2000. We are currently working on a paper to describe this algorithm as applied to viscous mixed element type grids but have not

published it because improvements are under development in regards to computational efficiency.