Turbulent Free-Surface Flows: A Comparison Between Numerical Simulations and Experimental Measurements
D.Dommermuth,1 M.Gharib,2 H.Huang,2 G.Innis,1 P.Maheo,2 E.Novikov,3 J.Talcott,1 D.Wyatt1 (1Science Applications International Corporation, 2California Institute of Technology, 3University of California at San Diego, USA)
The results of large-eddy simulations (LES) are compared to measurements of a model-scale frigate wake. The measurements are compared to two types of large-eddy simulations: a nonlinear free-surface approach and a low Froude-number approximation. An unique procedure has been developed for initializing three-dimensional LES using two-dimensional measurements of the mean and fluctuating velocity fields of a cross section of the wake. The procedure also works well for extending Reynolds-averaged Navier-Stokes simulations of near-field flows into the far wake. The numerical results show that turbulent roughening of the ocean surface scales according to three-dimensional similarity theory. An extensive inertial range is established in the kinetic and potential energy spectras. There is no evidence of the effects of two-dimensional turbulence in the free-surface elevation. Based on analyses of energy spectra, third-order moments, and probability distributions, the performance of the Smagorinsky model is satisfactory in the inviscid limit. Preliminary results of some recent grid-stirred turbulence and splitter-plate experiments are discussed. The numerical predictions of the evolution of the wake agree well with towing-tank measurements.
The results of towing-tank experiments are used to initialize and to validate large-eddy simulations (LES) of ship wakes. The towing-tank experiments were performed at the Maritime Research Institute, Netherlands (MARIN). The experiments were performed on a frigate. The model was propelled, and it was equipped with a rudder. Details of the experimental procedure are provided in Hoeskstra (1991).
Reed, et al (1990) and Hoeskstra (1991) review the physical attributes of ship wakes. In the turbulent portion of the wake, two pairs of vortices tend to form due to flow separation at the stern. A pair of bilge vortices forms due to the tapering of the parallel middle-body into the stern. A second pair of vortices forms near the free-surface due to the stern rising up to make room for the propeller. The bilge vortices can be entrained by propeller slipstreams and swirls. There can also be complex interactions with rudders. The free-surface vortices rise up into the drag wake of the hull and run parallel to the surface due to image vortex effects.
Both pairs of vortices are observed in the frigate experiments that had been performed at MARIN. The MARIN measurements are compared to two types of large-eddy simulations. The first type uses nonlinear free-surface boundary conditions and the second type uses a low Froude-number approximation to model free-surface roughness. The subgrid-scale (SGS) turbulence model uses the Smagorinsky model. Forcing at low wavenumbers is used to balance the dissipation at high wavenumbers to provide the proper dissipation rate based on experimental measurements. The free-surface code uses an adjustment scheme to prevent the generation of spurious high-frequency waves (see Dommermuth, 1994b). The free-surface code also implements a free-surface boundary-layer approximation that is valid for clean free surfaces (see Dommermuth, 1994a).
The LES codes are initialized using the experimental data. On a transverse cut all three velocity components are measured. The mean axial and transverse velocity components are used to assign the mean velocity in a cross section of fluid that is modeled using LES. LES is used to model the turbulent portion of the wake. The Kelvin wave disturbance is subtracted out of the measured data using linear wave-cut analysis. We superimpose on top of the mean velocity field a fluctuating velocity field that is based on the rms velocity fluctuations that are measured in the towing tank. The fluctuating velocity field is initialized using a realization of fully-developed turbulence with a –5/3 inertial range. We assume that the temporal behavior predicted in our numerical results is related to the spatial behavior measured in the laboratory experiments through a Galilean transformation. The numerical predictions of the evolution of the wake correlate well with towing-tank experiments.
Consider the unsteady incompressible flow of a Newtonian fluid under a free surface, and let ui(x,t)=(u,v,w) represent the three-dimensional velocity field as a function of time. Applying Helmholtz's decomposition gives
where (x,y,z,t) is a velocity-potential which describes the irrotational flow and Ui=(U,V,W) is a solenoidal-field which describes the vortical flow such that
where ∇2 is the Laplace operator. Since satisfies Laplace's equation and the divergence of the rotational field U is chosen zero, the total velocity field u conserves mass. The potential models the effects of waves, whereas the solenoidal velocity field U models the effects of shear currents and turbulence.
Based on this Helmholtz decomposition of the velocity field, the total-pressure Π is defined in terms of a vortical pressure P and an irrotational pressure, as follows:
Here, the velocity and pressure terms are respectively normalized by uc and where uc is the characteristic velocity and ρ is the density. is the Froude number and ℓc is the characteristic length. The vertical coordinate z is positive upward, and the origin is located at the mean free surface. Substituting these decompositions (1 & 4) into the Navier-Stokes equations gives
where Re=ucℓc/v is the Reynolds number and v is the kinematic viscosity. Mij is a subgrid-scale model and Si is a stirring force. The stirring term, which is a linear operator, is required to impose the correct dissipation rate for LES of wakes. The stirring coefficient α is calculated to provide a prescribed turbulent dissipation rate :
where V is the volume of fluid.
The divergence of the momentum equations (5) used in combination with the mass-conservation equations (2 & 3) can be used to derive a Poisson equation for the vortical pressure:
As shown by Dommermuth (1993), the vortical pressure is also subject to a solvability condition due to the imposition of Neumann boundary conditions.
Free-Surface Boundary Conditions
The Helmholtz decomposition of the velocity field requires that an additional boundary condition be imposed on the free surface. An expedient boundary condition that can be specified is that the normal component of the rotational velocity is zero on the free surface:
where z=η(x,y,t) is the free-surface elevation and ni is the unit-normal on the free surface. The preceding constraint imposed on the rotational velocity field means that the evolution of the free-surface elevation is entirely prescribed in terms of the free-surface elevation itself and the velocity potential, as follows:
where everything is evaluated on the exact position of the free surface, z=η.
The normal stress on the free surface must balance with the atmospheric pressure and the surface tension:
where Pa is the atmospheric pressure due to wind forcing, is the Weber number, and σ is the surface tension. d/dt=∂/∂t+ηt∂/∂z is a substantial derivative. In addition to the normal-stress condition there are also two tangential-stress conditions that are provided in Dommermuth (1993).
The normal stress condition (10) and the corresponding tangential stress conditions are appropriate only if the free-surface boundary-layer is sufficiently resolved. Since it is difficult to resolve the boundary-layer at full-scale Reynolds numbers, a boundary-layer approximation is required.
Longuet-Higgins (1991), Lundgren (1989), and Dommermuth (1994a) provide examples of free-surface boundary-layer formulations. For a clean free surface, the boundary-layer approximation of the normal stress condition is
This equation correctly predicts the laminar dissipation for small amplitude waves. The two
tangential stress conditions are not required in this boundary-layer approximation. Excluding the boundary-layer correction, the same free-surface boundary conditions are also used for the Euler equations. This approximation is valid when the temporal and spatial scales of the turbulence are longer than the scales of the waves.
The boundary-layer approximation models the dominate effects of refraction, radiation, scattering, and dissipation of waves by shear currents and turbulence through the combination of the kinematic condition for the rotational velocity (Equation 8) and the vortical pressure term (P) in the normal stress condition (Equation 11). The generation of vorticity by spatial and temporal changes in the free-surface elevation is not modeled, but the boundary-layer approximation does model the dissipation of waves due to the direct action of viscosity.
If the Froude number is sufficiently low, Dommermuth (1994b) shows that the free-surface elevation is hydrostatically balanced with the portion of the pressure that is induced by the vortical flow:
For a clean free surface, Dommermuth, et al (1994) show that the corresponding boundary conditions for the solenoidal velocity field are as follows:
These boundary conditions are satisfied on the plane z=0. At this level of approximation, the velocity potential, which represents dispersive wave behavior, is zero.
Based on the nonlinear free-surface conditions (8–11) and the low-Froude number approximation (12 & 13), two LES approaches have been developed. Aside from the treatment of the free surface, both numerical schemes are the same.
Sub-Grid Scale Modeling
The performance of several subgrid scale (SGS) models was discussed by Dommermuth & Novikov (1993). As pointed out by them, a major advantage of the Helmholtz decomposition is that SGS models for the turbulent portion of the flow can be developed independently of the potential portion of the flow. As a result, waves will not be excessively attenuated, and the potential flow will not adversely affect the turbulence model. A Smagorinsky SGS model is used here:
where Cs is the Smagorinsky coefficient, Ci is the coefficient of the isotropic turbulence model, Δ is the grid scale, δij is the Kronecker delta function, and Sij is the rate-of-strain tensor:
|S|=(SijSij)1/2 is the magnitude of the strain tensor. Most LES absorb the isotropic SGS stresses into the pressure-gradient term in the Navier-Stokes equations, but in free-surface flows this term is required to predict turbulent roughening of the free surface (see Dommermuth & Novikov, 1993).
Germano, et al's (1991) test filter approach in combination with Lilly's (1991) least-squares procedure is used to predict the Smagorinsky coefficient and the coefficient of the isotropic turbulence model. Details of the procedure are provided in Dommermuth & Novikov (1993).
The SGS terms associated with wavy portion of the velocity field and the cross-coupling between the wavy and vortical portions of the velocity field are not modeled. At low Froude numbers, the effects of these terms are not expected to be large because wave generation is minimal.
A stirring force is applied to maintain the proper dissipation rate (6). The stirring term forces the lowest wavenumbers, whereupon the energy cascades to higher wavenumbers. In wavenumber space, the stirring operator () is given by:
(κ)=((20+30 cos(κ)+ 12 cos(2κ)+2 cos(3κ))/64)4, (16)
where κ is a nondimensional wavenumber based on grid spacing and 0≤κ≤π. This filter is applied to the velocities along the three coordinate axes. For example, let x, y, and z represent low-wavenumber filtering along the x, y, and z axes, respectively. Let denote the three-dimensional Fourier transform of the Solenoidal velocity field, and let denote the Fourier transform of the stirring force, then
Similar filters are described in Dommermuth & Novikov (1993).
A Galilean approximation is used to relate the spatial evolution of a ship wake to the temporal evolution of a LES:
where x is the longitudinal distance downstream in the wake and t is the corresponding time in the LES. Uo is the ship speed.
Based on this Galilean approximation, we further assume that
where F is a physical quantity, an overbar denotes time averaging, brackets denote spatial averaging, L is the length of the LES, and T is the duration of a large-scale event in the wake. Xo and To are positions in space and time where the ship wake and LES correspond.
For LES of ship wakes, we assume that the potential and solenoidal velocity fields can be decomposed as follows:
where w is the wave disturbance due to Kelvin waves, u is the wave disturbance due to turbulence, ui is a mean disturbance due to turbulence, and vi is a fluctuating disturbance due to turbulence. We assume that the Kelvin waves do not interact strongly with the turbulent wake.
If the mean temporal velocities and the Reynolds stresses are provided on a transverse cut of the wake, the assumptions (18–20) can be used to initialize a LES. A cross section of the wake is modeled using LES. The first step in the procedure is to use linearized wave-cut analysis to represent the Kelvin wave disturbance in the wake:
where x and y are respectively the longitudinal and transverse positions in the wake. The angular dependence Φ(θ) is calculated using Neumann-Kelvin theory or wave-cut data. The water-particle velocity due to the Kelvin waves is calculated using Uw=∇w.
LES is used to simulate the vortical portion of the flow and the free-surface waves that are induced by the turbulence. The disturbance due to Kelvin waves is subtracted from the total velocity field. The remaining disturbance is primarily due to the effects of vortical flow. The vortical portion of the flow can be initialized using the results of RANS, experiments, or empirical data.
In the case of experimental data, the mean temporal disturbance due to the vortical flow is projected onto a solenoidal field using a
vorticity-stream function formulation. Let Ū denote the mean vector velocity with respect to time on a transverse cut of a wake, and let <u> denote the mean velocity with respect to the streamwise direction in the LES. In general, Ū is not solenoidal, nor does Ū satisfy the proper boundary conditions. To impose solenoidality and to enforce boundary conditions, we first calculate the vorticity (Ω):
This vorticity is used to calculate a vector stream function (Ψ):
where free-slip boundary conditions are initially imposed on the free surface and the bottom, and periodic boundary conditions are imposed in the horizontal plane. The mean velocity in the LES is calculated by definition:
which is solenoidal and satisfies the boundary conditions.
To complete the initial conditions, a fluctuating velocity field (vi) is superimposed over the mean velocity field in the LES. Let i denote a fully-developed realization of homogeneous anisotropic turbulence with zero mean, <i>=0. We project i onto vi using a deformation tensor:
where Bi is a solenoidal deformation tensor and Cd is a deformation constant.
Let denote the fluctuating velocity field as modeled by RANS or measured during experiments. Let denote the turbulent kinetic energy on a transverse plane of the wake, and let represent the LES counterpart. These two turbulent kinetic energies are fitted together using a least-squares and averaging procedure to determine the deformation constant:
where the surface integral is over a transverse cut of the wake S.
A natural choice for Bi, the deformation tensor, is the mean velocity field:
This choice of Bi is solenoidal, as required, and it leads to a fluctuating velocity field that satisfies the free-surface boundary conditions in the low-Froude number limit. Since this choice of Bi assigns a relationship between the large-scale features of the flow and the small-scale fluctuations, it could also potentially be used in RANS to provide closures for the Reynolds stresses.
Another choice for Bi is the measured rms velocity distribution:
By construction, this Bi is also solenoidal. If the condition ∂B1/∂z=0 on z=0 is imposed, then the free-surface boundary conditions will be satisfied by the fluctuating velocity field in the low Froude-number limit. The advantage of this approach is that it correlates well with the experimental input.
In summary, the mean and rms velocities on a transverse cut of a wake are used to seed a three-dimensional LES. The initial data can be provided by either experiments or RANS. A Galilean transformation is used to relate the duration of the simulation to the distance downstream in the experiments. Spatial averaging in the LES corresponds to temporal averaging in the experiments. First, the Kelvin wave disturbance is subtracted from the experiments, and then the mean velocities are projected onto a solenoidal field. Finally, a realization of fully-developed turbulence is fitted to the measured turbulent kinetic energy distribution on a transverse cut of the wake to construct the three-dimensional fluctuating velocity field in the LES.
The Navier-Stokes equations, and the boundary and initial conditions are discretized using 2nd-order finite differences. A fully-staggered grid is used in the numerical simulations. The momentum equations (5), kinematic condition (9), and the normal stress condition (11) are integrated with respect to time using a third-order Runge-Kutta scheme. The convective terms in the momentum equations (5) are expressed in conservation-law form. Each stage of the Runge-Kutta scheme is formulated to inhibit the accumulation of errors in the divergence of the rotational flow field. The rotational pressure is used to project the rotational velocity onto a solenoidal field (3 & 7) with zero normal velocity on the free surface (8). Laplace's equation for the potential (2) and Poisson's equation for the rotational pressure (7) are solved at each stage of the Runge-Kutta scheme, and a solvability condition is enforced for the rotational pressure. A multigrid solution scheme is used to solve the three-dimensional elliptic equations (see Dommermuth, 1994a). The z-coordinate is mapped onto a flat plane (see Dommermuth, 1994b). This mapping is applied to Laplace's equation (2), the Poisson equation for the pressure (7), the momentum equations (5), and all the boundary conditions. Periodic boundary conditions are used on the sides of the domain and free-slip boundary conditions are used on the bottom. For the low Froude-number approximation, free-slip boundary conditions are also used on the top (13). The numerical algorithms have been implemented on CM-5 parallel computers. The results of several validation studies are provided in Dommermuth & Novikov (1993).
Grid-stirred turbulence experiments are being performed in an open water tunnel with a width of 45.7cm and a water depth of 41cm. The free-stream velocity is Uc=28.36cm/s. Turbulence is generated by a mesh grid. The wire diameter is Do=1.5mm and the spacing between wires is Dw=8.3mm. Based on mesh spacing, the Reynolds number is Rem= UcDw/v=2350.
The wake of a vertical flat plate is being studied in a splitter-plate facility. The cross-section of the flow channel is 101.6cm wide and 53.65cm deep. The splitter plate is Dp= 7.6cm wide in its thickest section, and it tapers to zero at its end. The speed of the stream is Uc=40cm/s. Based on plate thickness, the Reynolds number is Rep=UcDp/v=30,400.
Velocity measurements are performed using Digital Particle Image Velocimetry (DPIV). DPIV provides instantaneous 2-D velocity fields in a plane (see Willert & Gharib, 1991). For the grid-stirred turbulence experiments, three planes have been measured: one vertical plane in the center of the water tunnel and two horizontal planes, 3mm and 110mm below the free surface. The area of each image plane is about 180 mm×130mm and the leading edge is about 235mm downstream of the wire mesh. 1000 images, which corresponds to 66.7 seconds of measurement, has been acquired in each plane to enable statistical analyses. Analysis of the rms velocity shows that averaging over 1000 samples gives relative errors less than one percent.
Towing-tank experiments were performed at MARIN using a model-scale frigate. The frigate was propelled, and it was equipped with a rudder (see Hoeskstra, 1991). Based on the towing speed (Uc=2m/s) and the length of the beam (ℓc=0.758m), which are used
to normalize the LES results, the Reynolds, Froude, and Weber numbers are respectively Re=1.52×106, and We= 4.10×104.
Data was measured at five transverse planes using laser doppler velocimetry. Longitudinal and transverse wave cuts were also performed. The velocity measurements were performed at about 200 points on each transverse plane. The sampling density was about 6 points across the draft and about 16 points across the beam. The sampling density was highest near the free surface and behind the stern and became coarser near the edges of the wake.
The experimental data has been interpolated onto a computational domain that is four beams wide and one beam deep. The length of the computational domain is one beam. Details of the two-dimensional interpolation procedure, the Kelvin-wave removal, and the noise reduction are provided in Innis (1993).
A fully-developed realization of homogeneous anisotropic turbulence is required to initialize LES of ship wakes. The resulting dataset is projected onto a transverse cut of the wake as discussed in Section 2.5. Once the LES has been initialized, it is advanced in time to simulate the spatial development of a ship wake. In the next two sections, we discuss LES of grid-stirred turbulence and ship wakes. The grid-stirred turbulence results, which are discussed in the first section, are used to initialize a ship wake study, which is discussed in the second section.
Numerical simulations of grid-stirred turbulence have been performed using the low Froude-number approximation. Figure 1 shows energy and surface-pressure spectra. The surface-pressure is the pressure on the plane z=0, where the free-surface approximation is made. Results are shown for 1283 and 2563 grid points. Stirring at low wavenumbers is used to maintain a constant kinetic energy with a rms velocity equal to one. The sides of the computational domain are one unit long. The kinematic viscosity is zero for these simulations. In this infinite Reynolds number limit, the only sink of energy is provided by the SGS model.
An energy cascade is established as indicated by the –5/3 and –7/3 slopes in the inertial ranges of the velocity and the surface-pressure spectra, respectively. The Kolmogorov constant for these simulations (Ko ≈ 2) agrees with other LES and DNS numerical simulations (see Chasnov 's review, 1991, Vincent & Meneguzzi, 1991), but it is higher than most experimentally observed values, which typically are about 1.5. Since the surface pressure can be directly related to the surface elevation through the relation provided in Equation (12), the –7/3 power-law behavior in the pressure spectra shows that the effects of two-dimensional turbulence are not likely to be observed in the far wake of a ship where the turbulence is well-stirred and the Froude number is low.
is fitted using a least-squares technique. For the 1283 simulation the fitting is performed for 20≤κ≤400, and for the 2563 simulation, 20≤κ≤800.
Table 1. Power-law behavior of energy spectra.
The results of other researchers seem to suggest that a free surface is capable of supporting a reverse energy cascade, which is a property of two-dimensional turbulence (see Sarpkaya's review, 1995). Observations of two-dimensional turbulence may be based on low-Reynolds-number flows that are not well-stirred. If the flow is well-stirred and the Reynolds number is sufficiently high, the present numerical results
indicate that the free surface does not support a reverse energy cascade.
Figure 2 shows two energy spectra measured in the grid-stirred turbulence experiments. One spectrum is based on data measured close to the free surface (z=–3mm) and the other is based on data measured in the bulk of the flow (z=–110mm). The data has been normalized by the free-stream velocity and the mesh spacing. For wavenumbers greater than the wavenumber at which turbulence is injected, ki>2.5, both spectra are close to a Kolmogorov –5/3 power-law. The presence of the free surface does not appear to influence the shape of spectra, which is consistent with the LES results. Due to limitations of the spatial resolution of DPIV in the present study, the energy spectra at very small scales, k>25 have not yet been resolved. A DPIV study that is zoomed in on the flow is currently underway.
Higher-order moments associated with the velocity increments are also useful measures of how well a turbulence model is performing. Following Novikov (1992), the velocity increment is defined as
where x is a vector position, r is a vector offset, and x′=x+r. The radial component of the velocity increment δui(r,x) is given by
where ni=rir–1, ri are the components of the offset, and r is the magnitude of the offset. This leads to Kolmogorov's formula for third-order moment, which is valid in the inertial range:
where the brackets denote volume averaging and is the dissipation rate.
Figure 3 shows the third-order moment normalized by its expected behavior in the inertial range. The normalized third-order moment asymptotically approaches one for small radial distances. The curve should be one throughout the inertial range. The Smagorinsky turbulence model appears to perform well for local disturbances, but its performance for long-range disturbances requires improvement. We note that the establishment of an extensive inertial range is adequate for remote-sensing applications, which only require good prediction of second-order moments. Applications include wave roughening, wave scattering, and turbulent wave dissipation. In order to predict intermittent turbulence effects, third and higher-order moments are required.
Figure 4 shows the probability distributions of the velocity increments for three radial off-sets. The velocity increments have been normalized such that their rms value is one. A Gaussian distribution is also included along with results from the field and laboratory experiments of Praskovsky & Oncley, 1994. The probability distributions are shown for various radial offsets. The radial offsets for the experiments are normalized by the Kolmogorov length-scale (ηk). The radial offsets in the LES are not normalized by ηk because dissipation scales are not readily available in LES.
The experimental results are strongly non-Gaussian, much more so than the numerics. However, the experiments do show a tendency to become more Gaussian as the radial offset increases. We note that the probability distribution can never become completely Gaussian because the third-order moment is non-zero in the inertial range. LES, with zero kinematic viscosity, corresponds to the limit r/ηk → ∞. In this context, the LES results represent a limiting value to the experiments, which is supported somewhat by the results illustrated in Figure 4 as r/ηk increases. The DNS results of Vincent & Meneguzzi (1991) show similar behavior. In the dissipation range, their DNS results are strongly non-Gaussian, but when the radial offset is outside the dissipation range and inside the inertial range, their results are more Gaussian.
A Model-Scale Frigate
The data from the MARIN experiments (see configuration 81, Hoekstra, 1991) measured at x=4.62 beams astern is used to initialize a LES using the procedure described in Section 2.5. The duration t=3.95 of the LES is long enough to advance the data to the next set of measurements at x=8.57. Based on the difference in kinetic energies between the two stations, a constant dissipation rate is assigned. This dissipation rate is maintained using stirring as described in Section 2.4. A more complex temporal behavior for the dissipation rate could be assigned, perhaps based on self-similarity arguments, but since the simulation is short, a constant dissipation rate is sufficient.
Table 2 provides the correlation coefficients for two different cases and several physical quantities. The two cases correspond to the two techniques for initializing the fluctuating portion of the velocity field. Case 1 uses the mean velocity field (see Equation 27), and Case 2 uses the rms velocity field (see Equation 28). In both cases, a low Froude-number approximation is used to model the wake.
Table 2. Correlation Coefficients.
Initially, the correlation coefficients are very high for all physical quantities, especially for the mean axial velocity (<U1>) and the mean axial vorticity (<Ω1>). For the initial conditions, these quantities just undergo a minor correction to ensure that the initial velocity field is solenoidal and to enforce the boundary conditions.
The correlation coefficient for the turbulent kinetic energy is also very high initially. The correlation coefficient for Case 1 is initially lower than for Case 2 for the turbulent kinetic energy. The mean velocity field that is used in Case 1 tends to emphasize regions where the magnitude of the velocity is high. However, most turbulence production occurs in regions where the gradients of the mean velocity field are high. There are gradients of the mean velocity in the fitting procedure (see Equation 25), but these terms tend to be dominated by the terms containing gradients of the grid-stirred turbulence field. Case 2, which uses the rms velocity as a fit, tends to do better in the high-gradient regions. These issues are more evident in the color Figures 5 and 6, which will be discussed later. A combination of mean velocity and mean vorticity may be used to provide a better fit to the turbulent kinetic energy field.
At the next station, at time t=3.95, the correlation coefficients for the mean axial velocity and the turbulent kinetic energy are high, whereas the coefficient for the mean axial vorticity is low. The low correlation coefficient in the vorticity field is due to the breakup of a pair of bilge vortices.
Case 1 has also been run using the nonlinear free-surface code. Comparisons to the low Froude-number approximation indicate that wave generation is not very strong. The correlation coefficients between the hydrostatic approximation (see Equation 12) and the nonlinear free-surface elevation vary from 0.8 near the middle of the simulation to 0.5 at the end of the simulation. The low-Froude approximation appears to be reasonable four beams astern, and further back in the wake, where the turbulence is weaker, it would get even better.
We had earlier hypothesized that when the mean velocity is used as a deformation tensor, Equations 25 and 27 could be used in RANS to provide closures for the Reynolds stresses. To test this hypothesis, we fitted a fluctuating velocity field to the stations at x=4.62 and x=8.57 in the MARIN experiments. The deformation constant (Cd) changed by less than 15% between the two stations. This provides some evidence that the initialization procedure could also be used as a RANS closure.
measurements and the numerical predictions for the two cases. Figure 5 is at t=0 and Figure 6 is at t=3.95. The mean turbulent kinetic energy, mean streamwise velocity, and mean axial vorticity are respectively shown in the (a) top, (b) middle, and (c) bottom plots.
For the mean axial velocity, the drag wake of the hull (blue) is located near the surface, and the thrust wake of the propeller (red) is concentrated below the free surface. Some spreading of the wake is observed between times t=0 and t=3.95. In both cases, LES predicts the evolution of the axial velocity very well.
The plots of mean axial vorticity show that the free-surface vortices are concentrated in a thin layer near the free surface. They spread laterally outward. The pair of bilge vortices have been entrained in the slipstream of the propeller and split by the rudder. The experimental data is ragged at t=0 because there are only a few sample points across the wake. Between time t=0 and t=3.95, the bilge vortices break up and the free-surface vortices grow weaker. The free-surface vorticity also spreads out laterally. The spreading of the wake is more evident in the free-surface elevation, which is not shown. The two numerical cases appear to perform equally well.
As illustrated in the plots of turbulent kinetic energy (tke), high concentrations of tke are initially located in the region between the drag and thrust wakes of the experiments. The numerics do not seem to capture this high concentration very well. As discussed earlier, it may partially be attributed to the gradients of the fluctuating velocity field dominating the gradients of the mean quantities that are used in the fitting procedure. The effect of using different realizations of turbulence to initialize the flow should also be investigated.
A new procedure for initializing LES and simulating wakes has been developed. The procedure provides a framework for testing SGS models. Based on this procedure, two types of LES capabilities have been developed: a nonlinear free-surface capability and a low Froude-number approximation. These capabilities have been used to simulate grid-stirred turbulence and the wake of a model-scale frigate.
Neither the results of LES or preliminary grid-stirred turbulence experiments shows evidence of a reverse energy cascade. For well-stirred flows at sufficiently high Reynolds numbers, the numerically-predicted and experimentally-measured energy spectra scale according to three-dimensional similarity theory, both in the bulk of the flow and near the free surface. The roughening of the free surface also scales according to three-dimensional similarity theory.
Comparisons between LES predictions and towing-tank measurements are very good. As predicted by LES, the attenuation of the mean axial-velocity field and the turbulent kinetic energy, and the breakup of the axial vorticity field agrees with measurements. The enhanced spreading of the wake near the free surface is also predicted well by LES. Additional validation studies are currently underway using the results of grid-stirred turbulence and splitter-plate experiments.
DGD, GEI, JCT, and DCW are supported by ONR under contract number N00014–93-C-0046. MG, HTH, and PM are supported by ONR under contract number ONR-N00014–92-J-1610. EAN is supported under contract numbers ONR-N00014–92-J-1610 and ONR-14–94– 1-0040. Dr. Edwin P.Rood is the program manager. This work is also supported in part by the Army Research Office contract number DAALO3–89-C0038 with the University of Minnesota Army High Performance Computing Research Center (AHPCRC) and the DoD Shared Resource Center at the AHPCRC. The numerical simulations have been performed on the CM-5 computers at Naval Research Laboratory and the AHPCRC.
 Chasnov, J.R. ( 1991) Simulation of the Kolmogorov inertial subrange using an improved subgrid model. Phys. of Fluids, 3(1), 188–200.
 Dommermuth, D.G. ( 1993) The Laminar Interactions of a Pair of Vortex Tubes with a Free Surface . J. of Fluid Mech., 246, 91– 115.
 Dommermuth, D.G. ( 1994a) Efficient simulation of short- and long-wave interactions with applications to capillary waves. J. Fluids Eng., 116, 77–82.
 Dommermuth, D.G. ( 1994b) The initialization of vortical free-surface flows. J. Fluids Eng., 116, 95–102.
 Dommermuth, D.G., Novikov, E.A., & Mui, R.C.Y. ( 1994) The interaction of surface waves with turbulence. In Proc. of the ASME Symp. on Free-Surface Turbulence, Lake Tahoe.
 Dommermuth, D.G. & Novikov, E.A. ( 1993) Direct-numerical and large-eddy simulations of turbulent free-surface flows. In Sixth Inter. Conf. on Numer. Ship Hydro., Iowa City, 239–270.
 Germano, M., Piomelli, U., Moin, P., & Cabot, W.H. ( 1991) A dynamic subgrid-scale eddy viscosity model. Phys. Fluids, A 3(7), 1760–1765.
 Hoekstra, M. ( 1991) Macro wake features of a range of ships. Maritime Research Institute, Netherlands, MARIN Report No. 410461–1-PV.
 Innis, G.E. ( 1993) Computations of full-scale wakes from model-scale data. Science Applications International Corporation, SAIC Report No. SAIC-93/1139.
 Lilly, D.K. ( 1991) A proposed modification of the Germano subgrid-scale closure method . Phys. Fluids, A 4(3), 633–635.
 Longuet-Higgins, M.S. ( 1991) Theory of weakly damped Stokes waves: a new formulation and its physical interpretation. J. Fluid Mech., preprint.
 Lundgren, T.S. ( 1989) A free-surface vortex method with weak viscous effects. In Mathematical Aspects of Vortex Dynamics, ed. by R.E.Caflisch, SIAM, 68–79.
 Monin, A.S. & Yaglom, A.M. ( 1975) Statistical Fluid Mechanics II. The MIT Press, Cambridge, MA
 Mui, R.C.Y. & Dommermuth, D.G. ( 1995) The vortical structure of a near-breaking gravity-capillary wave. J. Fluids Eng., 117, 355–361.
 Novikov, E.A. ( 1992) Probability distribution for three-dimensional vectors of velocity increments in turbulent flow. Phys. Rev. A, 46(10), 6147–6149.
 Praskovsky, A. & Oncley, S. ( 1994) Probability density distribution of velocity differences at very high Reynolds numbers. Phys. Rev. Lett., 73(25), 3399–3402.
 Reed, A.M., Beck, R.F., Griffin, O.M., & Peltzer, R.D. ( 1990) Hydrodynamics of remotely sensed surface ship wakes. In SNAME Trans., 98, 319–363.
 Sarpkaya, T. ( 1995) Vorticity, free surface, and surfactants. Ann. Rev. of Fluid Mech., 28, 83–128.
 Smagorinsky, J. ( 1963) General circulation experiments with the primitive equations. Mon. Weather Rev., 91, 91–164.
 Vincent, A. & Meneguzzi, M. ( 1991) The spatial structure and statistical properties of homogeneous turbulence . J. Fluid Mech., 225, 1–20.
 Willert, C. & Gharib, M. ( 1991) Digital particle image velocimetry. Exp. in Fluids, 10, 181–183.