Parameterization of the Near-Field Internal Wave Field Generated by a Submarine
James W. Rottman, Kyle A. Brucker, Douglas Dommermuth, Dave Broutman
aa r X i v : . [ phy s i c s . f l u - dyn ] O c t th Symposium on Naval HydrodynamicsPasadena, California, 12-17 September 2010
Parameterization of the Near-Field Internal Wave FieldGenerated by a Submarine
James W. Rottman , Kyle A. Brucker , Douglas Dommermuth and DaveBroutman Science Applications International Corporation, USA Computational Physics, Inc., USA
ABSTRACT
We attempt to gain some insight into the modeling ofthe generation of internal waves produced by submarinestraveling in the littoral regions of the ocean with the useof high fidelity numerical simulations. These numericalsimulations are shown to be capable of simulating highReynolds number flow around bodies, including the ef-fects of stable stratification. In addition, we use the re-sults of these detailed numerical studies to test and revisethe source distribution parameterizations of the near-fieldwaves that have been used in analytical studies based onlinear theory. Such parameterizations have been shownto be useful in initializing ray-tracing schemes that canbe used for computing wave propagation through real-istic oceans with variable background properties. Forsimplicity, we focus on the idealized case of a spheri-cal body traveling horizontally at constant speed througha uniformly stratified fluid.
INTRODUCTION
The motion of a submarine through a stratified oceanproduces an internal wave field that could be used, ifaccurately forecasted, for detection and tracking. Thesewaves are produced by the vertical displacement of thefluid as it flows over the submarine body and by the dis-turbance of the fluid by the motions in the submarine’swake. The wake motions consist of turbulent eddies andthe bulk motion due to the collapse of the partially mixedwake region towards its equilibrium density level. In lit-toral regions of the ocean, where stratification is strongand submarines travel slowly, the body-generated wavespredominate. In the open ocean, where stratification gen-erally is weaker, the wake-turbulence generated wavesare dominant.A very fast and accurate model has beendeveloped to compute the propagation of internalwaves through realistic ocean environments. Themodel, a modification of the methods described in Broutman and Rottman (2004), Rottman et al. (2006)and Rottman and Broutman (2008) involves a ray andcaustic solution in the Fourier transform domain, whichis mapped into a spatial solution by inverse Fourier trans-form. This is a more practical approach than calculatingthe ray and caustic solution directly in the spatial do-main and is general enough to treat background flowswith depth dependent shear and stratification.The propagation model requires initial condi-tions for the rays. These initial conditions can be pro-vided by a near-field approximation of the generationof the waves. A convenient approximation is to repre-sent the waves emitted by a horizontally moving bodyas that due to a horizontally moving and vertically oscil-lating distribution of sources in a depth-dependent back-ground. The horizontally moving distribution of sourcesmodels the body-generated waves and the vertical oscil-lation of the source is an approximation of the unsteadywave generation by the turbulent wake. This idealizationof the generation of the internal wave field was proposedby Voisin (1994) and Dupont and Voisin (1996) and isbased on some ideas from the laboratory experiments ofGorostov and Teodorovich (1981, 1982, 1983). Similarideas have been used by Robey (1997) for an ocean witha particularly sharp thermocline. However, these theoriesfor the parameterization of the near-field internal wavefield, and particularly the wave field generated by the tur-bulent wake, have not been thoroughly validated.The focus of this paper is to use the results ofnumerical simulations to test and revise the source dis-tribution parameterizations of the near-field waves pro-duced by submarines traveling in the littoral regions ofthe oceans. The study is restricted to the idealized caseof a sphere traveling horizontally at constant speed U through a vertically stratified fluid at low Froude number, Fr , defined as F r = U/ ( N a ) in which a is a measure ofthe radius of the body and N is the (constant) buoyancyfrequency of the fluid.The numerical simulations are for the stratifiedflow over and around the body so that the internal waveeld due to the displacement of the fluid by the body isknown. The simulations are fully nonlinear, and use atechnique similar to the Cartesian-grid free-surface cap-turing code (NFA) of Dommermuth et al. (2007). Thenear wake solution is coupled with the technique ofDommermuth et al. (2002), which efficiently simulatesturbulent wakes to distances in excess of 1000 body di-ameters down stream, using the temporal approximationto simulate the internal waves generated by the collapseof the turbulent wake. Both the body generated as wellas the wake generated waves can be simulated. This al-lows a direct comparison between the numerically sim-ulated wave field and the near-field internal wave fieldcomputed from the parameterized source distribution.The Fourier-ray model, as specialized for thecase of uniform stratification, is described in the nextsection, including the oscillating source distribution rep-resentation for flow produced by the body. The numer-ical model is outlined in following section, along witha number of validation studies. The comparison of theFourier-ray model results with the numerical simulationsis described in the final section. THE FOURIER-RAY MODEL
The main aspects of the Fourier-ray model for hori-zontally homogeneous background conditions that mayhave vertical variations of the buoyancy frequencyand currents are described in Rottman et al. (2004b),Rottman et al. (2004a), Broutman et al. (2003) andBroutman and Rottman (2004). Here our focus is on theparamterization of the generation of linear internal wavesby the flow of a stratified fluid over and in the wake ofa sphere. For this purpose, we consider the simpler caseof linear internal waves generated by a body moving hor-izontally at a steady speed U through a stably stratifiedBoussinesq fluid with constant buoyancy frequency andno background currents. As a means of parameterizingthe waves generated by a turbulent wake, we will oscil-late the sphere vertically at a constant frequency σ , asdiscussed in more detail later in this section.The coordinate system is r = ( x, y, z ) , with z positive upwards, and is fixed to the mean position of thebody. In this reference frame, the background flow is U = ( U, , . The background buoyancy frequency isa constant value N . We solve for the vertical displace-ment η ( r , t ) . Any other dependent variable can be com-puted from the vertical displacement using standard lin-ear wave polarization relations, as shown in Gill (1982).An equation for η ( r , t ) can be derived from thelinearized Boussinesq equations of motion for a uniformbackground, as described for example by Miles (1971),Lighthill (1978), or Broutman and Rottman (2004), [ D ∂ z + ( D + N )( ∂ x + ∂ y )] η = D∂ z q. (1) in which D = ∂ t + U ∂ x . The source distribution rep-resenting the body is specified in the form q ( r , t ) = q ( r ) e − iσt , in which σ is the oscillation frequency of thesource distribution.The solution to this equation can be found byFourier transform in the horizontal directions. η ( r , t ) = e − iσt ∞ Z Z −∞ ˜ η ( k, l, z ) e i ( kx + ly ) dk dl . (2)where ˜ η ( k, l, z ) is the vertical eigenfunction, ˜ η = sgn( z )2 πi ˆ ωmQ e imz /B m . (3)We are using subscript notation for partial derivatives, e.g. B m = ∂B/∂m , and m is the internal wave verticalwavenumber given by the linear dispersion relation, m = ± ( k + l ) / ( N / ˆ ω − / . (4)in which the internal wave wavenumber k = ( k, l, m ) and the intrinsic frequency ˆ ω is given by ˆ ω = σ − kU. (5)The function B is, see Lighthill (1978), B = ˆ ω m − ( N − ˆ ω )( k + l ) . (6)Finally, Q ( k ) is the three-dimensional Fouriertransform of q . Note that in Q , m is treated as a func-tion of k, l through the dispersion relation.As the factor e − iσt accounts for all of the time-dependence in the present model, this solution can beconsidered as the long-time limit of an initial value prob-lem in which the motion is started from rest and the bodyasymptotically in time oscillates vertically with constantfrequency σ .R EFLECTIONS FROM UPPER AND LOWERBOUNDARIES
To account for wave reflections from upper andlower horizontal boundaries, if they exist, we followthe method of Broutman et al. (2003), which we outlinehere. For simplicity, we restrict attention to a backgroundthat is depth independent. Each time a ray returns to anyfixed depth z above the body in a channel of total depth H after reflecting once from the top of the channel andonce from the bottom, the wave phase mz has changedby an amount m H . To account for the effects of thisand subsequent reflections, we multiply ˜ η ( k, l, z ) in (3)by the sum S = ∞ X n =0 e in mH (7) = ie − imH / mH ) . (8)2his is a divergent sum, since the individual terms do notvanish as n → ∞ , however the sum can be evaluatedin the sense of generalized functions (see Eq. (1.2.2)of Hardy (1949). In any case, the sum diverges when mH = π , when the interference of the reflected waves isperfectly constructive. We eliminate this divergence byadding a small damping factor in the form of an imag-inary wavenumber for the vertical wavenumber, or bylimiting S to a finite number of terms. The number ofterms can be chosen to represent the number of reflec-tions at a given time of interest, for given k, l , as deter-mined by a group velocity calculation. We have experi-mented with both methods but have used only the secondmethod for the results presented here.S OURCE DISTRIBUTION
We consider internal waves generated by thevertically oscillating sphere of radius a in a uniformflow of speed U . In the limit of large Froude number F r = U/N a , as shown by Gorostov and Teodorovich(1982) and Dupont and Voisin (1996), the flow associ-ated with this motion can be represented by the followingsource distribution function q ( r , t ) = − h U xa + hσe − iσt za i δ ( r − a ) , (9)where δ is the Dirac delta function. The vertical dis-placement amplitude of the oscillation of the sphere is h ,and its vertical velocity is hσe − iσt . Its three-dimensionalFourier transform Q is Q ( k , t ) = 34 iπ − a (cid:2) U k + hσe − iσt m (cid:3) j ( Ka ) Ka , (10)where K = | k | , and j ( z ) = ( sin z ) /z − ( cos z ) /z is thespherical Bessel function of order unity.The non-oscillating portion of the source dis-tribution is an accurate representation of non-stratifiedflow over a sphere, but has been shown (and we willshow here) that it also is a very good representation ofthe flow over sphere even for moderate to low Froudenumbers. The oscillating portion of the source dis-tribution has been used to model internal wave gener-ation by turbulent eddies in the wake of an obstacle(Dupont and Voisin (1996) based on some experimen-tal work by Gilreath and Brandt (1985)). Following theguidelines of Dupont and Voisin (1996), the dominantinternal waves generated by the eddies in the turbulentwake are simulated by choosing a value of . for theStrouhal number St = aσ/piU and a source frequencyof σ = N πStF r . The amplitude of the vertical oscilla-tion of the source, represented by the factor h in (9), is atthis stage chosen to match the amplitude of the observedwaves. N UMERICAL PROCEDURE
Putting this altogether, the spatial solution isfound at any specified height z by evaluating the ray so-lution (3), with B given by (6) and Q given by (10), andmultiplying the result by (8), and finally taking the in-verse Fourier transform. THE NUMERICAL MODEL
The computer code Numerical Flow Analysis (NFA),Dommermuth et al. (2007), originally designed to pro-vide turnkey capabilities to simulate the free-surface flowaround ships, has been extended to have the ability toperform high fidelity stratified sub-surface calculations.The governing equations are formulated on a Cartesiangrid thereby eliminating complications associated withbody-fitted grids. The sole geometric input into NFA isa surface panelization of the body. No additional grid-ding beyond what is used already in potential-flow meth-ods and hydrostatics calculations is required. The easeof input in combination with a flow solver that is im-plemented using parallel-computing methods permit therapid turn around of numerical simulations of high- Re stratified fluid-structure interactions.The grid is stretched along the Cartesian axesusing one-dimensional elliptic equations to improve res-olution near the body. Away from the body, where theflow is less complicated, the mesh is coarser. Details ofthe grid stretching algorithm, which uses weight func-tions that are specified in physical space, are provided inKnupp and Steinberg (1993).G OVERNING E QUATIONS
Consider a turbulent flow in a stratified fluid.Physical quantities are normalized by characteristic ve-locity ( U o ), length ( L o ), time ( L o /U o ), density ( ∆ ρ ), andpressure ( ρ o U o ) scales, where ρ o is the reference densityand ∆ ρ is the change in density over the characteristiclength scale. Let ρ and u i respectively denote the nor-malized density and three-dimensional velocity field asa function of normalized space ( x i ) and normalized time( t ). The conservation of mass is ∂ρ∂t + ∂u j ρ∂x j = 0 . (11)For incompressible flow with no diffusion, ∂ρ∂t + u j ∂ρ∂x j = 0 . (12)Subtracting (12) from (11) gives a solenoidal conditionfor the velocity: ∂u i ∂x i = 0 . (13)3or an infinite Reynolds number, viscousstresses are negligible, and the conservation of momen-tum is ∂ρu i ∂t + ∂∂x j ( ρu j u i ) = − ρ o ∆ ρ (cid:20) ∂p∂x i + ρRi B δ i + τ i (cid:21) , (14)where p is the normalized pressure and τ i is a normal-ized stress that will act tangential to the surface of thebody. δ ij is the Kronecker delta function. The sub-grid scale stresses are modeled implicitly O’Shea et al.(2008). Ri B is the bulk Richardson number defined as: Ri B ≡ ∆ ρρ o gL o U o , (15)where g is the acceleration of gravity. The bulk Richard-son number is the ratio of buoyant to inertial forces.The normalized density is decomposed in termsof the constant reference density plus a small departurewhich is further split into a mean and a fluctuation: ρ = ρ o ∆ ρ + ρ ( x ) + ˜ ρ ( x i , t ) . (16)Here, ρ ( x ) and ˜ ρ ( x i , t ) are respectively the normalizeddensity stratification and the normalized density pertur-bation. The pressure, p is then decomposed into the dy-namic, p d , and hydrostatic, p h , components as: p = p d + p h . (17)The hydrostatic pressure is defined in terms of the refer-ence density and the density stratification as follows: ∂p h ∂x i = − ( ρ o ∆ ρ + ρ ) Ri B δ i . (18)The substitution of (16) into (11) provides anew expression for the conservation of mass: ∂ ˜ ρ∂t + ∂∂x j ( u j ( ρ + ˜ ρ )) = 0 . (19)The substitution of (16)-(18) into (14) and us-ing (11) to simplify terms gives a new expression for theconservation of momentum: ∂u i ∂t + ∂∂x j ( u j u i ) = −
11 + γ ( ρ + ˜ ρ ) (cid:20) ∂p d ∂x i + ˜ ρRi B δ i + τ i (cid:21) , (20) where γ = ∆ ρ/ρ is the characteristic density differencedivided by the reference density. If γ << , a Boussi-nesq approximation may be employed in the precedingequation to yield: ∂u i ∂t + ∂∂x j ( u j u i ) = − ∂p d ∂x i − ˜ ρRi B δ i − τ i . (21)If the background stratification is linear, we let ∆ ρ = − L o dρ d dx d , (22)where a supersript d denotes a dimensional variable. Fora linear background stratification, Ri B = F r − o , where F r o = U o ( N L o ) − is the internal Froude number and N is the Brunt-V¨ais¨al¨a frequency defined as: N = − gρ o dρ d dx d , (23)The bulk Froude number, a ratio of inertial to gravita-tional forces, is defined as: F r b ≡ U o ( gL o ) / . (24)Multiplying Ri B and F r b , (15) times (24) squared,yields the ratio of buoyant to gravitational forces: Ri B F r o = ∆ ρρ = γ. (25)When the background stratification is linear, γ = F r b /F r o .The momentum equations using either (20) or(21) and the mass conservation equation (19) are inte-grated with respect to time. The divergence of the mo-mentum equations in combination with the solenoidalcondition (13) provides a Poisson equation for the dy-namic pressure. The dynamic pressure is used to projectthe velocity onto a solenoidal field and to impose a no-flux condition on the surface of the body. The details ofthe time integration, the pressure projection, and the for-mulation of the body boundary conditions, are describedin the next three sections.T IME I NTEGRATION
A second-order Runge-Kutta scheme is used tointegrate with respect to time the field equations for thevelocity and density. During the first stage of the Runge-Kutta algorithm, a Poisson equation for the pressure issolved: ∂∂x i (cid:20)
11 + γ ( ρ + ˜ ρ k ) (cid:21) ∂p kd ∂x i = ∂∂x i (cid:18) u ki ∆ t + R ki (cid:19) , (26)4here R ki denotes the nonlinear convective, buoyancy,and stress terms in the momentum equation, (20), at timestep k . u ki , ˜ ρ k , and p kd are respectively the velocity com-ponents, density, and dynamic pressure at time step k . ∆ t is the time step. For the next step, this pressure isused to project the velocity onto a solenoidal field. Thefirst prediction for the velocity field ( u ∗ i ) is u ∗ i = u ki + ∆ t (cid:18) R ki − (cid:20)
11 + γ ( ρ + ˜ ρ k ) (cid:21) ∂p kd ∂x i (cid:19) . (27)The density is advanced using the mass conservationequation (19): ˜ ρ ∗ = ˜ ρ k − ∆ t ∂∂x j (cid:2) u kj ( ρ + ˜ ρ k ) (cid:3) (28)A Poisson equation for the pressure is solved again dur-ing the second stage of the Runge-Kutta algorithm: ∂∂x i (cid:20)
11 + γ ( ρ + ˜ ρ ∗ ) (cid:21) ∂p ∗ d ∂x i = ∂∂x i (cid:18) u ∗ i + u ki ∆ t + R ∗ i (cid:19) , (29) u i is advanced to the next step to complete one cycle ofthe Runge-Kutta algorithm: u k +1 i = 12 (cid:0) u ∗ i + u ki (cid:1) + ∆ t (cid:18) R ∗ i − (cid:20)
11 + γ ( ρ + ˜ ρ ∗ ) (cid:21) ∂p ∗ d ∂x i (cid:19) , (30)and the density is advanced to complete the algorithm: ˜ ρ k +1 = ˜ ρ ∗ + ˜ ρ k − ∆ t ∂∂x j (cid:2) u ∗ j ( ρ + ˜ ρ ∗ ) (cid:3) . (31)E NFORCEMENT OF N O -F LUX B OUNDARY C ONDITIONS
A no-flux condition is satisfied on the surface ofthe submerged body using a finite-volume technique. u i n i = v i n i , (32)where n i denotes the unit normal to the body that pointsinto the into the fluid and v i is the velocity of the body.Cells near the body may have an irregular shape, depend-ing on how the surface of the body cuts the cell. Let S b denote the portion of the cell whose surface is on thebody, and let S o denote the other bounding surfaces ofthe cell that are not on the body. Gauss’s theorem is ap-plied to the volume integral of (26): Z S o + S b ds (cid:20) n i γ ( ρ + ˜ ρ k ) (cid:21) ∂p kd ∂x i = Z S o + S b ds (cid:18) u ki n i ∆ t + R i n i (cid:19) . (33) Here, n i denotes the components of the unit normal onthe surfaces that bound the cell. Based on (27), a Neu-mann condition is derived for the pressure on S b as fol-lows: (cid:20) n i γ ( ρ + ˜ ρ k ) (cid:21) ∂p kd ∂x i = − u ∗ i n i ∆ t + u ki n i ∆ t + R i n i . (34)The Neumann condition for the velocity (32) is substi-tuted into the preceding equation to complete the Neu-mann condition for the pressure on S b : (cid:20) n i γ ( ρ + ˜ ρ k ) (cid:21) ∂p kd ∂x i = − v ∗ i n i ∆ t + u ki n i ∆ t + R i n i . (35)This Neumann condition for the pressure is substitutedinto the integral formulation in 33: Z S o ds (cid:20) n i γ ( ρ + ˜ ρ k ) (cid:21) ∂p kd ∂x i = Z S o ds (cid:18) u ki n i ∆ t + R i n i (cid:19) + Z S b ds v ∗ i n i ∆ t . (36)This equation is solved using the method of fractional ar-eas. Details associated with the calculation of the areafractions are provided in Sussman and Dommermuth(2001) along with additional references. Cells whose cutvolume is less than 2% of the full volume of the cell aremerged with neighbors. The merging occurs along thedirection of the normal to the body. This improves theconditioning of the Poisson equation for the pressure. Asa result, the stability of the projection operator for the ve-locity is also improved (see Equations 27 and 30).E NFORCEMENT OF N O -S LIP B OUNDARY C ONDITIONS
The stress τ i is used to impose partial-slip andno-slip conditions on the surface of the body using abody-force formulation as follows: τ i = β ( u i − v i ) δ ( x − x + b ) , (37)where β is a body-force coefficient, v i is the velocityof the body, x is a point in the fluid, and x + b is a pointslightly outside the body. Equation (37) forces the fluidvelocity to match the velocity of the body. Note that free-slip boundary conditions are recovered with β = 0 andno-slip boundary conditions are imposed as β → ∞ .Dommermuth et al. (1998) discuss modelingusing body-force formulations. Recently, there have5een several studies which use a similar body boundaryconditions in both finite volume (Meyer et al. (2010a),Meyer et al. (2010b) and finite element (Hoffman(2006a), Hoffman (2006b), John (2002) simulations.The finite element simulations of Hoffman (2006b) atvery high Reynolds numbers are able to predict the liftand drag to within a few percent of the consensus valuesfrom experiments, using a very economical number ofgrid points. The finite volume implementations have alsoshown promise albeit at more modest Reynolds numbers( Re = 3 , ).In the following sections it is shown that finitevolume simulations with Re → ∞ are able to accu-rately predict flow about bluff-bodies using a partial-sliptype boundary condition to model the effects of the un-resolved turbulent boundary layer. Validation
As discussed in O’Shea et al. (2008) the advective termsin the momentum equations are handled with the fluxbased limited QUICK scheme of Leonard (1997). Here,this treatment is extended to include the advective termsin the density equation. In this section three types ofnumerical simulations are performed to assess the valid-ity of different aspects of the computational model pro-posed in the preceding section. First, a canonical turbu-lent shear layer is simulated to test the ability of NFA toaccurately resolve highly turbulent flows, while remain-ing stable in the absence of explicit turbulence models, ormolecular viscosity. Second, the unstratified flow overa sphere is considered to test the ability of NFA to ac-curately predict flow separation on a curved bluff body.Lastly, the stratified flow over a sphere is considered atseveral values of internal Froude number ranging from . − . to ensure that the effects of density stratifica-tion are accurately captured.T URBULENT S HEAR L AYER
The first type of simulation used to validateNFA is a temporally evolving shear layer. An example ofthis type of flow is shown in Figure 1 to provide a visualreference to the amount of disorder these flows generate,and a reference to the problem setup.A shear layer develops when two streams offluid with different velocities are brought together. Thistype of flow is representative of many engineering andgeophysical flows. The analysis of this type of flowis simplified because the statistical description of theflow reduces to a single dimension. It is often consid-ered a canonical flow to study fundamental aspects ofturbulence, and has been extensively studied both exper-imentally (Corcos and Sherman (1976), Bell and Mehta
Figure 1:
NFA simulation of a turbulent mixing layer. Theillustration is an x − z plane of instantaneous streamwise, x ,velocity. Note the wide range in the scales of motion, which istypical of shear generated turbulence. (1990) and Spencer and Jones (1971), amongst oth-ers) and numerically (Rogers and Moser (1994),Pantano and Sarkar (2002) and Brucker and Sarkar(2007), again amongst others). These detailed studiesprovide a rich database to compare with, and thesecomparisons serve to assess the ability of NFA tosimulate highly turbulent flows.Here, the velocity streams, U and U , are ofequal magnitude and opposite sign, with the velocitydifference across the two streams being ∆ U ≡ U − U . The model flow evolves temporally in a domainthat moves with the mean convective velocity U C =( U + U ) / . The thickness of the mixed regiongrows in time and can be characterized by the vorticitythickness, δ ω ( t ) and then momentum thickness, δ θ ( t ) ,respectively defined as: δ ω ( t ) ≡ ∆ U ( d h u i /dx ) max , (38) δ θ ( t ) = ∞ Z −∞ (cid:18) − h u i / ∆ U (cid:19) dx . (39)In the preceding equations the h·i operator denotes anaverage over the x − x plane, and fluctuations withrespect to the mean are obtained by applying a Reynoldsdecomposition, u i = h u i i + u ′ i . Initial Conditions
The initial conditions are constructed as a meancomponent plus a small random disturbance. The ad-justment procedure of Dommermuth et al. (2002), inwhich the mean is held constant until the productionreaches a steady value, is used. The adjustment timewas δ ω, / ∆ U . Figure 2 shows the h u ′ u ′ i cross-correlation responsible for turbulent production. By t =8 the a constant value was reached, and at t = 10 the6 x / δ ω,0 < u u > / ∆ U t= 0.331t= 1.584t= 3.169t= 7.628t= 7.076t= 7.352t= 7.907 Figure 2:
Reynolds shear stress responsible for turbulent pro-duction at several times during the adjustment phase. mean was allowed to evolve. This type of adjustmentprocedure removes the arbitrary addition of “turbulentlike” fluctuations, and allows for repeatable initial con-ditions. The mean initial conditions are: h u i = − ∆ U (cid:18) x δ ω, (cid:19) , h u i = h u i = h p i = 0 . (40)where δ ω, = 1 and ∆ U = 1 .The amplitude of the random fluctuations addedto the mean velocity was . U . The computationaldomain [ L , L , L ] = [36 . , . , . δ ω, was dis-critized with [ N , N , N ] = [384 , , computa-tional nodes. The time step was ∆ t = 0 . δ ω, / ∆ U .Periodic boundary conditions were used in the stream-wise ( x, x ) and span-wise ( y, x ) direction for all vari-ables. In the cross-stream direction ( z, x ), Dirichletboundary conditions were used for the vertical veloc-ity and Neumann boundary conditions were used for theother velocity components and pressure.S HEAR L AYER T URBULENCE
Figures 3(a)-(c) compare the Reynolds stresses, R ij = (cid:10) u ′ i u ′ j (cid:11) for i = j , obtained here with data fromprevious studies. The Reynolds stresses were evalu-ated by averaging the self-similar form at the followingtimes: δ ω, /U t = δ ω, / ∆ U . Thepeak streamwise, transverse, and spanwise turbulent in-tensities, shown in Figures 3(a)-(c), √ R / ∆ U = 0 . , √ R / ∆ U = 0 . , and √ R / ∆ U = 0 . agree wellwith previous DNS and experimental data. The experi-mental data shows a scatter of about with respect (a) -1 -0.5 0 0.5 1 x / δ ω (t) R / / ∆ U B&S 2007P&S 2002R&M 1994B&M 1990S&J 1971NFA (b) -1 -0.5 0 0.5 1 x / δ ω (t) R / / ∆ U B&S 2007P&S 2002R&M 1994B&M 1990S&J 1971NFA (c) -1 -0.5 0 0.5 1 x / δ ω (t) R / / ∆ U B&S 2007P&S 2002R&M 1994B&M 1990S&J 1971NFA
Figure 3: (a) Streamwise, (b) Spanwise, (c) Cross-stream r.m.s. velocities. The abbreviations in the leg-end are as follows: B&S 2007 refers to the DNS ofBrucker and Sarkar (2007), P&S 2002 refers to the DNS ofPantano and Sarkar (2002), R&M 1994 refers to the DNS ofRogers and Moser (1994), B&M 1990 refers to the experimentsof Bell and Mehta (1990), and S&J 1971 refers to the experi-ments of Spencer and Jones (1971). to the simulation data. The shape of the self-similar pro-files also agrees well. Comparisons of the the Reynoldsstresses, R ij = (cid:10) u ′ i u ′ j (cid:11) for i = j , and the production of7a) -1 -0.5 0 0.5 1 x / δ ω (t) R / / ∆ U B&S 2007P&S 2002R&M 1994B&M 1990S&J 1974NFA (b) -4 -2 0 2 4 x / δ θ (t) P ∆ U / δ θ ( t ) NFAB&S 2007P&S 2002R&M 1994
Figure 4: (a) Reynolds shear stress h u ′ u ′ i . (b) Turbulent pro-duction. References in legend as in Figure (3). turbulent kinetic energy defined as: P ≡ − (cid:10) u ′ i u ′ j (cid:11) ∂ h u i i ∂x j = − h u ′ u ′ i ∂ h u i ∂x , (41)in the current simulations to that from an incompress-ible DNS (Rogers and Moser (1994),Brucker and Sarkar(2007) and compressible low Mach number DNS(Pantano and Sarkar (2002) are shown in Figure 4(a)-(b),respectively. The results from the current simulationsagree well with previous DNS and laboratory data.In summary, the turbulent shear layer simula-tions have been successfully bench marked against pre-vious laboratory and DNS data, and show NFA’s abilityto correctly capture the energetically important scales ofturbulent motion without requiring explicit modeling. Flow over a Sphere
NFA simulations of flow over a sphere in uniform andstratified fluids are respectively discussed. The com-putational domain is shown schematically in Figure 5.
Figure 5:
A schematic of the computational domain.
The size of the computational domain [ L , L , L ] , thenumber of computational nodes [ N , N , N ] , the nearbody grid spacing, ∆ x NBi , the time step, ∆ t , the inter-nal Froude number, F r , and the near-body parameter β are provided in Table 1 for all simulations subsequentlydiscussed.All of the simulations used In-flow and Out-flow boundary conditions in the stream-wise ( x ) direc-tion, and Neumann boundary conditions in the span-wise( x ) direction. In the cross-stream direction ( x ), Dirich-let boundary conditions are used for the vertical velocity, u , and Neumann boundary conditions are used for theother velocity components, pressure, and density.In the NFA simulations the flow is acceleratedfrom rest by multiplying the free-stream current with thefollowing function: f ( t ) = n − exp h − ( t/T cur ) io . (42) T cur = 0 . was used in all simulations.U NIFORM FLUID
Contours of the instantaneous streamwise ve-locity on the x − z plane located at y = 0 are shownin Figure 6 for cases F r ∞ with β = 1 (top) and F r ∞ with β = 0 (bottom). When β = 0 the boundary con-ditions on the body are no-slip and no flow separationoccurs. When β = 1 the flow separates just shy
120 deg ,the value reported in Achenbach (1972). Here, a com-parison of the drag coefficient to theory and experimentsis used to assess the quality of the flow near the body.Since in the unstratified case the drag is determined bythe fields in the immediate vicinity of the body it is anexcellent metric for determing the quality of those fields.The drag coefficient, is defined as: C D ≡ F D / ρ U A P , (43)8 igure 6: NFA simulations of unstratified flow over a sphere. The illustrations are of the instantaneous stream-wise velocity on the x − z plane located at y = 0 . Top case F r ∞ with β = 1 . Bottom case F r ∞ with β = 0 . β . All normalized with the body diameter D and free-stream velocity U . The internal Froude number, F r , for each case is given by the case name. Cases labeled
F rX denote a domain used for runs at different Froudenumbers. Case L x L y L z N x N y N z ∆ t β ∆ x NBi
F r
19 20 10 512 384 384 0.0025 0.67 0.03
F r
19 20 5 512 384 256 0.0020 0.33 0.03
F r
19 20 10 1024 768 768 0.00125 1.00 0.01
F r
19 20 5 512 384 256 0.002 0.33 0.03
F r ∞
10 8 8 1024 768 768 0.001 1.00 0.008
F r ∞
10 8 8 1024 768 768 0.001 0.00 0.008
F rX C
12 6 6 384 256 256 0.005 0.67 0.03
F rX M
12 12 12 384 384 384 0.005 0.67 0.03
Time C D Theory F = m du/dtAchenbach (1972)Fr ∞ (β =0.00)Fr ∞ (β =1.00)Startup Quasi-Steady Figure 7:
Drag coefficient for flow over a sphere at high-Reynolds number. The references in the legend are: Achenbach(1972) refers to the experiments of Achenbach (1972) while thetheoretical prediction refers to equation 44. where F D is the drag force and A P is the project frontalarea of the body.The drag force on the body is calculated by inte-grating the normal pressure over the surface. Achenbach(1972) reports that the viscous contribution to the drag is ∼ of the total at Re = 5 × , and hence the the totaldrag should be well approximated by the pressure drag athigh- Re . The viscous contribution is not calculated.The coefficient of drag for cases F r ∞ with β = 1 and F r ∞ with β = 0 are shown in Figure 7. Asnoted in the discussion of Figure 6 the free-slip bound-ary conditions, used in case F r ∞ , did not allow for flowseparation and hence there was almost complete pressurerecovery in the lee of the sphere and hence no drag. Forcase F r ∞ , in which separation occured at the correctlocation, the drag coefficient was . the same value re- ported in the experiments of Achenbach (1972).During the startup phase of the simulations, f ( t ) is non-zero and there is an added mass componentof the drag: F AMD = mU dfdt (44)The drag due to the added mass is also shownin Figure 7, and the agreement between both simulationsand the prediction is excellent.F LOW OVER A SPHERE IN A STRATIFIED FLUID
Having validated the drag on a sphere in thecase without stratification, attention is now turned tothe case with stratification. To proceed comparisonsare made between existing theoretical and experimen-tal data sets. The stratification of the fluid can causebody generated waves and significant drag at low Froudenumbers ( ≤ . ), while at moderate Froude numbers( . ≥ F r ≤ ) the potential energy gained over-coming the poles of the sphere aids in pressure recov-ery in the lee and a small drag reduction. Following,Lofquist and Purtell (1984) the change in drag due tostratification, ∆ C D , is defined as: ∆ C D ( F r ) = C D ( F r ) − C D ( F r ∞ ) . (45)Figure 8 shows the change in the drag coeffi-cient, ∆ C D , as a function of the internal Froude num-ber, F r = U/ (2 rN ) . NFA simulations at F r =0 . , . , . , . , . , . , . , . , . are comparedwith the theories of Greenslade (2000), Voisin (2007),and Gorodtsov and Teodorovich (1982) along with theexperiments of Lofquist and Purtell (1984), Vosper et al.(1999), Shishkina (1996), and Mason (1977).In summary, Figure 8 shows that NFA is ableto predict correctly the change in drag for both lowand moderate Froude numbers, and appears to correctly10 .5 1 1.5 2 Fr ∆ C D FrX C FrX M Vosper et al (1999) HS2Shishkina (1996) Mason et al (1977)L&P N=0.205L&P N=0.324L&P N=0.457L&P N=0.644L&P N=0.787L&P N=0.891L&P N=1.091L&P (1984)G&T (1982)Voisin (2007)Grnsld (2000) C=3.0Grnsld (2000) C=2.5
Figure 8:
Change in the drag coefficient as a function of inter-nal Froude number. In the legend,
F rX C and F rX M denotecases with the coarse and medium resolutions reported in Ta-ble 1. predict the transition from low to high Froude numberregimes, which occurs between F r ≈ . and F r ≈ . ,although this prediction currently is based on only onecomputed data point. INTERNAL WAVES
We now compare the internal waves computed by NFAfor stratified flow over a sphere with the field predictedby linear wave theory, using an oscillating source distri-bution to parameterize the generation of the waves bythe sphere and its wake, as discussed in The Fourier-Ray Model section. As this is a preliminary study, weconsider flows for just two moderate Froude numbers,
F r = 1 and . The sphere is towed in the positive x di-rection, and each plot shown below is either a horizontal ( x, y ) cross section at a height of two sphere diameters( z = 2 D ) above the centerplane of the sphere or a ver-tical ( x, z ) plane through the wake centerline. The plotsshow the x component u of the wavefield velocity.For the Fourier-ray calculations, the inverse Figure 9:
Along track vertical ( x, z ) plane at y = 0 of the nu-merical solutions at the wake center for F r = 1 (upper panel)and
F r = 4 (lower panel). The plotted variable is u/U , the x -component of the velocity of the wavefield. Fourier transform (2) was approximated discretely on awavenumber grid of 1024 by 512 in k, l , for
F r = 1 and
F r = 4 . In addition to resolving the flow features, theextent of the wavenumber grid must be chosen to limitperiodic wrap-around errors, which result from the dis-crete approximation of the inverse Fourier transform. Fora single depth, the theoretical solution takes only a fewseconds to calculate on a standard PC.Figure 9 shows the horizontal velocity u/U on avertical plane through the wake centerline ( y = 0 ) for thetwo cases F r = 1 (top pane) and
F r = 4 (bottom pane).The lower Froude number case shows a well-defined leewave field as well as a weakly turbulent wake. The leewave field is steady with respect to the sphere, but down-stream of about x/D = − internal waves generated bythe turbulent wake can be seen. For the higher Froudenumber, the lee waves are very weak and the turbulentwake is stronger and shows evidence of periodic fluctua-tions. Figure 10 show results for the case of F r =1 . The numerical simulation result is shown in the top11anel, and the linearized theory, with σ = 0 , in thelower panel. In this case the upper and lower bound-aries of the simulations are sufficiently far away for thereto be no reflections that would be seen in the domainshown in these figures. The two results are very close,with some small differences appearing downstream of x/D = − near y = 0 . At this relatively low Froudenumber the body-generated waves dominate over thewake-generated waves. It appears from these results thatthe proposed source representation of the sphere pro-duces an accurate internal wave field, even though thesource function we have used, (9), is most accurate forthe representation of a solid sphere only in the limit ofhigh Froude number. However, the discrepancies that areseen further downstream are due to weak internal wavesgenerated by the wake turbulence that can be seen in fig-ure 9. This effect is better visualized in figure 11, whichis a plot of u/U in the vertical plane located at y = 0 .Only the positive z half of the plane is shown. In theseplots the linear theory (lower panel) shows waves ema-nating only from the sphere, where as the numerical sim-ulations clearly show waves emanating from the sphereand from locations downstream of the sphere.Figure 12 show results for the case of F r = 1 ,now with the upper and lower boundaries at z = ± D ,which are close enough for the reflected waves to affectthe observed wavefield. The numerical simulation resultis shown in the top panel, and the linearized theory, with σ = 0 , in the lower panel. Again, the two results arevery close, with some small differences appearing down-stream of x/D = − near y = 0 , although the numeri-cal simulations show some reflections from the sidewallboundaries. The linear theory was computed withoutsidewall boundaries. The discrepancies seen downstreamand near the centerline are attributable to wake-generatedwaves, which are not accounted for in the linear theory.Figures 13 and 14 show u/U on a horizontalplane and a vertical plane, respectively, for F r = 4 . Herethe linear theory, with σ = 0 , is a poor representation ofthe solution, though it does give a reasonable predictionfor the lateral extent of the wavefield. This is the Froudenumber regime in which wake-generated internal wavesare expected to dominate over the body-generated waves.The wave amplitudes shown are quite small, and appar-ently significantly affected by the wake-generated waves.These results are very preliminary; we intend to continuethe study by extending the numerical simulations fartherdownstream and comparing these results with the lineartheory that includes combinations of nonzero source os-cillation frequencies. The hope is that we will be able todetermine the correct distribution of oscillation frequen-cies and oscillation amplitudes so that the linear theorywill accurately reproduce the internal wave fields in thesesimulations. x / D y / D −15 −10 −5 0−505 −0.04−0.0200.020.040.06x / D y / D −15 −10 −5 0−505 −0.04−0.0200.020.040.06 Figure 10:
A comparison for
F r = 1 from the numerical sim-ulations (upper panel) and from linear theory (lower panel), onthe ( x, y ) plane located at z = 2 D . The plotted variable is u/U , the x -component of the velocity of the wavefield. CONCLUSIONS
In this paper, we have extended the capabilities of thenumerical model NFA to high Reynolds number flowsaround obstacles in a stratified fluid. We have demon-strated that NFA (1) is capable of accurately reproducingthe physics associated with highly turbulent flows; (2) iscapable of accurately reproducing highly turbulent flowswhich separate due to the geometry of the body obstruct-ing the flow; and (3) is capable of simulating stratifiedflows at low and high Froude numbers.This work is another step in seeing how wellray theory can simulate the internal wavefield generatedby stratified flow past an obstacle. In a previous paper,Rottman et al. (2004a) we compared ray theory with lab-oratory experimental results. In this paper we have be-gun a comparison of the ray simulation with the numer-ical model results for uniform background with reflect-12 / D z / D −15 −10 −5 02468 −0.1−0.0500.05x / D z / D −15 −10 −5 02468 −0.1−0.0500.05 Figure 11:
A comparison for
F r = 1 from the numerical sim-ulations (upper panel) and from linear theory (lower panel), forthe positive − y region of ( x, y ) plane located at z = 2 D . Theplotted variable is u/U , the x -component of the velocity of thewavefield. ing upper and lower boundaries. The results show thatthe source distribution produces a very good representa-tion of the internal wave field for low Froude numbersnear unity where we would expect waves generated bythe body itself to dominate. At higher Froude numbers,where the dominated wave generation is by the turbu-lent wake, the comparisons are not as good. This studywill be continued in the future, now that NFA is availablefor these kinds of comparisons, to determine specificallywhat distribution of oscillating sources best models thewave generation by eddies in the wake. Acknowledgements
This research was sponsored by Dr. Ron Joslin at theOffice of Naval Research (contract number N00014-08-C-0508), Dr. Tom C. Fu at the Naval Surface Warfare
Figure 12:
A comparison for
F r = 1 , with upper and lowerhorizontal boundaries located at z = ± D , from the numericalsimulations (upper panel) and from linear theory (lower panel).The plotted variable is u/U , the x -component of the velocityof the wavefield. References
Achenbach, E., “Experiments on the flow past spheres atvery high reynolds number,” J. Fluid Mech., Vol. 54,no. 3, 1972, pp. 565–575.Bell, J. H. and Mehta, R. D., “Development of atwo-stream mixing layer from tripped and untripped13 / D y / D −15 −10 −5 0−505 −10−8−6−4−202x 10 −3 x / D y / D −15 −10 −5 0−505 −2−1.5−1−0.500.511.5x 10 −3 Figure 13:
A comparison for
F r = 4 from the numerical sim-ulations (upper panel) with the linear theory (lower panel), onthe ( x, y ) plane located at z = 2 D . The plotted variable is u/U on a horizontal plane at z/D = 2 . boundary layers.” AIAA J., Vol. 28, 1990, pp. 2034–2042.Broutman, D. and Rottman, J., “A simplified fouriermethod for computing the internal wavefield gener-ated by an oscillating source in a horizontally moving,depth-dependent background,” Phys. Fluids, Vol. 16,2004, pp. 3682–3689.Broutman, D., Rottman, J. W., and Eckermann, S. D.,“A simplified fourier method nonhydrostatic mountainwaves,” J. Atmos. Sci., Vol. 60, 2003, pp. 2686–2696.Brucker, K. A. and Sarkar, S., “Evolution of an initiallyturbulent stratified shear flow.” Phys. Fluids, Vol. 19,no. 10, 2007, p. 105,105.Corcos, G. M. and Sherman, F. S., “Vorticity concentra-tion and the dynamics of unstable free shear layers,” J.Fluid Mech., Vol. 73, 1976, pp. 241–264. x / D z / D −15 −10 −5 056789 −2−1012x 10 −3 x / D z / D −15 −10 −5 056789 −2−1012x 10 −3 Figure 14:
A comparison for
F r = 4 from the numerical sim-ulations (upper panel) with the linear theory (lower panel), forthe positive − z region of the ( x, z ) plane located at y = 0 . Theplotted variable is u/U . Dommermuth, D., Innis, G., Luth, T., Novikov, E.,Schlageter, E., and Talcott, J., “Numerical simulationof bow waves,” Proceedings of the 22nd Symposiumon Naval Hydrodynamics, Washington, D.C., 1998,pp. 508–521.Dommermuth, D. G., O’Shea, T. T., Wyatt, D. C., Rat-cliffe, T., Weymouth, G. D., Hendrickson, K. L., Yue,D. K. P., Sussman, M., Adams, P., and Valenciano,M., “An application of cartesian-grid and volume-of-fluid methods to numerical ship hydrodynamics,”Proceedings of the 9 th International Symposium onNumerical Ship Hydrodynamics, held 5 – 8 August2007 in Ann Arbor, Michigan, 2007.Dommermuth, D. G., Rottman, J. W., Innis, G. E., andNovikov, E. A., “Numerical simulation of the wake ofa towed sphere in a weakly stratified fluid,” J. FluidMech, Vol. 473, 2002, pp. 83–101.14upont, P. and Voisin, B., “Internal waves generatedby a translating and oscillating sphere,” Dyn. Atmos.Oceans, Vol. 23, 1996, pp. 289–298.Gill, A. E., Atmosphere-Ocean Dynamics, AcademicPress, 1982.Gilreath, H. E. and Brandt, A., “Experiments on the gen-eration of internal waves in a stratified fluid,” AIAAJ., Vol. 23, 1985, pp. 693–700.Gorodtsov, V. A. and Teodorovich, E. V., “Study of in-ternal waves in the case of rapid horizontal motionof cylinders and spheres,” Fluid Dynamics, Vol. 17,1982, pp. 893–898.Gorostov, V. A. and Teodorovich, E. V., “On the gen-eration of internal waves in the presence of uniformstraight-line motion of local and nonlocal sources,”Izv. Atmos. Oceanic Phys., Vol. 16, 1980, pp. 699–704.Gorostov, V. A. and Teodorovich, E. V., “Two-dimensional problem for internal waves generated bymoving singular sources,” Fluid Dynamics, Vol. 16,1981, pp. 219–224.Gorostov, V. A. and Teodorovich, E. V., “Study of in-ternal waves in the case of rapid horizontal motionof cylinders and spheres,” Fluid Dynamics, Vol. 17,1982, pp. 893–898.Gorostov, V. A. and Teodorovich, E. V., “Radiation ofinternal waves by periodically moving sources,” Appl.Mech. Tech. Phys., Vol. 24, 1983, pp. 521–526.Greenslade, M. D., “Drag on a sphere moving horizon-tally in a stratified fluid,” J. Fluid Mech., Vol. 418,2000, pp. 339–350.Hardy, G. H., Divergent Series, Clarendon, 1949.Hoffman, J., “Adaptive simulation of the subcritical flowpast a sphere.” J. Fluid Mech., Vol. 568, 2006a, pp.77–88.Hoffman, J., “Simulation of turbulent flow past bluffbodies on coarse meshes using general galerkinmethods: drag crisis and turbulent euler solutions.”Comput. Mech., Vol. 38, 2006b, pp. 390–402.John, V., “Slip with friction and penetration with resis-tance boundary conditions for the navier-stokes equa-tions - numerical tests and aspects of the implementa-tion,” J. Comp. Appl. Math., Vol. 147, 2002, pp. 287–300.Knupp, P. M. and Steinberg, S., Fundamentals of gridgeneration, CRC Press, 1993. Leonard, B., “Bounded higher-order upwind multi-dimensional finite-volume convection-diffusion al-gorithms,” W. Minkowycz and E. Sparrow, eds.,Advances in Numerical Heat Transfer, Taylor andFrancis, Washington, D.C., 1997, pp. 1–57.Lighthill, J., Waves in Fluids, Cambridge UniversityPress, 1978.Lofquist, K. E. B. and Purtell, P., “Drag on a sphere mov-ing horizontally through a stratified liquid,” J. FluidMech., Vol. 148, 1984, pp. 271–284.Mason, P. J., “Forces on spheres moving horizontally inrotating stratified fluid,” Astrophys. Fluid Dyn., Vol. 8,1977, pp. 137–154.Meyer, M., Devesa, A., Hickel, S., Hu, X. Y., andAdams, N. A., “A conservative immersed interfacemethod for large-eddy simulation of incompressibleflows,” J. Comput. Phys., Vol. in press.Meyer, M., Hickel, S., and Adams, N. A., “Assessment ofimplicit large-eddy simulation with a conservative im-mersed interface method for turbulent cylinder flow,”Int. J. Heat Fluid Flow, Vol. in press.Miles, J. W., “Internal waves generated by a horizontallymoving source,” Geophys. Astrophys. Fluid Dyn.,Vol. 2, 1971, pp. 63–87.O’Shea, T. T., Brucker, K. A., Dommermuth, D. G., andWyatt, D. C., “A numerical formulation for simulatingfree-surface hydrodynamics,” Proceedings of the 27thSymposium on Naval Hydrodynamics, Seoul, Korea,2008.Pantano, C. and Sarkar, S., “A study of compressibilityeffects in the high-speed turbulent shear layer usingdirect simulation,” J. Fluid Mech., Vol. 451, 2002, pp.329–371.Robey, H. F., “The generation of internal waves by atowed sphere and its wake in a thermocline,” Phys.Fluids, Vol. 9, 1997, pp. 3353–3367.Rogers, M. M. and Moser, R. D., “Direct simulation ofa self-similar turbulent mixing layer.” Phys. Fluids,Vol. 6, 1994, pp. 903–923.Rottman, J. W. and Broutman, D., “Practical imple-mentation of two-turning point theory for model-ing internal wave propagation in a realistic ocean,”Proceedings of the 27 th International Symposium onNaval Hydrodynamics, held 5 – 10 October 2008 inSeoul, Korea, 2008.Rottman, J. W., Broutman, D., Spedding, G., and Di-amessis, P., “A model for the internal wavefieldproduced by a submarine and its wake in the lit-toral ocean,” Proceedings of the 26 th International15ymposium on Naval Hydrodynamics, held 17 – 22September 2006 in Rome, Italy, 2006.Rottman, J. W., Broutman, D., Spedding, G., and Meu-nier, P., “The internal wave field generated by the bodyand wake of a horizontally moving sphere in a strati-fied fluid,” Proceedings of the 15 th Australasian FluidMechanics Conference, held 13 – 17 December 2004in Sydney, Australia, 2004a.Rottman, J. W., Broutman, D., Spedding, G., andMeunier, P., “Internal wave generation by a hor-izontally moving sphere at low Froude number,”Proceedings of the 25 thth