1 INTRODUCTION

Internal waves generated by a submerged moving body in a stratified fluid represent one of the significant branches of stratified flow dynamics, and are closely related to military applications. In general, the waves can be classified into three types. The first is generated by a combined volume flux effect contributed by a body and its attached flow structure, referred to as a lee wave. The second is generated by gravitational collapse of a mixed region in the body wake, referred to as a wake-collapse internal wave. The third is generated by the collapse of large-scale turbulent structures in the body wake, referred to as a wake random internal wave (Wei et al., 2009). A literature survey shows that the first and third types have received considerable attention but reports about the second type remain scarce (Gilreath and Brandt, 1985; Bonneton et al., 1993; Wei et al., 2009; Brandt and Rottier, 2015). Summarizing the available research results shows that laboratory studies only confirm the generation of the second type of internal waves, without detailed analysis of their properties (Schooley and Stewart, 1963; Gilreath and Brandt, 1985). Even numerical studies do not directly consider the essential effects of both body geometry and propeller forcing (Chernykh et al., 2004). Thus, it is still unclear whether wake-collapse internal waves can be generated by a propeller-driven body in oceanic stratified conditions, and what the wave properties are if they are generated. To answer these two questions, the present study simulates turbulent wakes generated by a propeller-driven body in a linearly stratified ocean, and explores the generation and properties of wake-collapse internal waves.

The paper is organized as follows. A newly developed numerical model and corresponding solution method are given in Section 2. Model initialization is described in Section 3. Section 4 contains numerical results with different body speeds and stratification magnitudes. Conclusions are given in Section 5.

2 NUMERICAL MODEL AND NUMERICAL SOLUTION METHOD 2.1 Numerical model

The newly developed numerical model consists of Reynolds-averaged Navier-Stokes equations, a propeller body-force model, and turbulence closure model.

2.1.1 Reynolds-averaged Navier-Stokes equations

The Reynolds-averaged Navier-Stokes equations are used to describe the situation in which a propellerdriven body moves with constant speed in an initially quiescent medium. A similar situation, in which a thin vertical barrier moves uniformly in such a medium, has been studied by Houcine et al. (2012). From their numerical results, they identified transient internal waves and stationary lee waves. In the present study, the body was geometrically similar to the classic Suboff model constructed by the Defense Advanced Research Projects Agency in the USA. The Suboff model is composed of an axisymmetric hull, a sail, and four stern appendages. To mimic a realistic submerged body, the size of the model is enlarged 25 times, to have an overall length of 108.9 m and maximum diameter of 12.7 m. In practice, propellerinduced flow is three-dimensional, highly unsteady and turbulent, so a very small time step and spatial resolution are required to capture the detailed flow field around a propeller. This results in a very high computational cost (Phillips et al., 2009). Such treatment is necessary when propeller performance and efficiency are the only interest. However, when the turbulent wake of a propeller-driven body is the interest, the propeller effect is normally represented by a body-force model (Schetz and Favin, 1979; Stern et al., 1988; Paterson et al., 2003; Alin et al., 2010). The common method of implementing a body-force model is directly adding source terms that represent the propeller forcing in the momentum equation. The relevant body-force model used herein is discussed in Section 2.1.2.

The momentum equation uses the Boussinesq approximation; that is, the effect of density variation is only considered in the gravity term. Further, the density of seawater is assumed to be a function of temperature only, and the Coriolis force is neglected. Then, the equations of conservation of momentum, conservation of mass, temperature transport (Afanasyev and Peltier, 2001; Vlasenko et al., 2005), and state (Millero and Poisson, 1981) are

Here, u, v, and w are the velocity components corresponding to x, y, and z coordinates; t is time; T is temperature; ρ is density; ρ_a is the initial stable density profile; ρ₀ is the reference density; g is gravitational acceleration; F_x, F_y, and F_z are the components of force generated by a propeller; μ_h and μ_z are horizontal and vertical eddy viscosity coefficients for momentum; and К_h and К_z are the horizontal and vertical eddy diffusivity coefficients for temperature.

2.1.2 Propeller body-force model

In the model, the effect of the propeller is represented by forces F_x, F_y, and F_z. According to the equation proposed by Paterson et al. (2003), these forces in Eqs.1.1-1.3 can be alternatively expressed as

where T_p represents the propeller's total thrust; Q_p represents its total torque; R_h is the hub radius; R_p is the propeller radius; L is the chord length at 0.7R_p; and and (x₀, y₀, z₀) are coordinates of the propeller center.

The determination of T_p and Q_p in Eqs.5-7 is as follows. It is known that T_p is equivalent to the total drag F that a body experiences when moving at constant velocity without the propulsion of a propeller. Therefore, T_p can be indirectly obtained by calculation of F, which is calculated by

where C_d is the total drag coefficient; U is the body speed; and S is the total surface area. The determination of C_d is via steady-state numerical simulation without considering propeller forcing in a homogenous fluid. Thus, substituting C_d into Eq.8 makes F known and thereby obtains T_p. After T_p is available, Q_p is obtained via the following relationship (Yang et al., 2013).

where K_t is the thrust coefficient and K_q the torque coefficient.

2.1.3 Turbulence closure model

The effects of eddy viscosity and diffusivity are included in the newly developed model. A constant coefficient turbulence model (Legg and Adcroft, 2003; Fringer et al., 2006; Chen et al., 2014) has been used in many numerical simulations of oceanic internal waves and is also used here because of its simplicity and efficiency. The characteristic horizontal eddy viscosity coefficient μ_h was 10^-1 m²/s in the light of numerical experiments in a homogeneous fluid with the same grid. Because the horizontal eddy viscosity and diffusivity coefficients mainly depend on grid resolution (Vlasenko and Hutter, 2002), we set them to μ_h=К_h=10^-1 m²/s. In addition, invoked by the turbulence model used by Vlasenko and Stashchuk (2007), we took the vertical eddy viscosity and diffusivity coefficients as μ_z=К_z=1.5×10^-2 m²/s.

2.1.4 Boundary conditions

To ensure that the generated internal waves could not reach the boundaries in the period of interest, the computational domain of the model was dimensioned as x∈(-200, 800) m, y∈(-150, 150) m, and z∈(-300, 0) m. The center of the submerged body was located at (100.2, 0, -100) m and its front and tail were at x=45.7 m and x=154.6 m, respectively (Fig. 1). To keep the propeller-driven body fixed in the computational domain, a moving coordinate system was adopted (Guan et al., 2012; Rojanaratanangkule et al., 2012).

In the experiments, no normal flux conditions were applied to the boundaries, except the one on the left. Inflow and outflow boundary conditions were specified at the left and right boundaries, respectively. Free slip boundary conditions were applied to the top and bottom boundaries, and no-slip boundary conditions at the body surface. At the side boundaries, normal velocity gradients were taken as zero. Numerical boundary conditions were applied to pressure at each boundary.

2.2 Numerical solution

The aforementioned numerical model was solved using the open-source computational fluid software OpenFOAM. This software uses the finite-volume method and has been widely accepted by investigators for studying complex problems (Lin and Song, 2012; Flore et al., 2013; Stringer et al., 2014). In our study, discretization of the computational domain was accomplished by the SnappyHexMesh utility within the software, and the CheckMesh utility was used to guarantee that the generated grid was adequate for accurate computation. For the discretization of equations, the temporal term was discretized by the first-order implicit Euler scheme. The advection terms were discretized by a second-order upwind differencing scheme and the diffusion terms by a second-order central differencing scheme. The discretized pressure and velocity were coupled with the pressure implicit with splitting of operators algorithm. To ensure the dominance of physical dispersion over numerical dispersion and fine spatial resolution to capture key processes, a uniform step of 5 m was used for the background mesh. Then, the body surface, propeller and wake regions were refined by factors of five, six and two, respectively. To make the computation satisfy the Courant-Friedrichs-Lewy stability condition, an adaptive time-stepping technique was used to maintain the maximum Courant number at no larger than 0.5 during the entire numerical calculation.

3 MODEL INITIALIZATION

A detailed description of the parameter setup and initialization of the temperature profile in the newly developed model are given in this section.

3.1 Parameters in propeller body-force model

To fit the enlarged body model, a new, improved propeller model was constructed based on standard seven-blade marine propeller INSEAN E1619 (Chase and Carrica, 2013), for which the primary parameters were R_p=3.275 m, R_h=0.68 m, and L=0.092 m. The location of the propeller center was (152.3, 0, -100) m, according to that prescribed in the Suboff model. We took K_t=0.2 and K_q=0.04 per the INSEAN E1619 open-water performance curve (Chase and Carrica, 2013), thereby relating T_p and Q_p by

where T_p depends on U. In general, the speeds of a submerged moving body range from 1.5 to 9 m/s. Thus, we set U to 1.5, 4.5, and 7.5 m/s, which represent low-speed, middle-speed and high-speed situations, respectively. Corresponding to each speed, the T_p and Q_p were obtained by steady-state numerical experiments without the propulsion of a propeller. Here, when U=1.5 m/s, T_p=15.08 m⁴/s² and Q_p=19.76 m⁵/s²; when U=4.5 m/s, T_p=119.94 m⁴/s² and Q_p=157.12 m⁵/s²; when U=7.5 m/s, T_p =319.46 m⁴/s² and Q_p=418.50 m⁵/s².

3.2 Initial temperature profile

The model temperature profile was given by

where H=300 m is total water depth; T₁ is temperature at the sea surface; and T₂ is temperature at the set bed. Considering that the order of oceanic buoyancy frequency is 10^-3/s to 10^-2/s, we designed three temperature profiles, which represent weakly, moderately and strongly stratified oceans, respectively. Buoyancy frequencies of the three cases at the location of the body z=-100 m were 3.27×10^-3/s, 8.66×10^-3/s and 1.38×10⁻²/s, respectively.

4 MODEL RESULT

The input parameters of the designed numerical experiments are listed in Table 1, including the body speed, stratification magnitude, and moving depth. We use Case1 to expound the generation process and basic properties of wake-collapse internal waves, and compare Case1, Case2 and Case3 to address the effect of body speed on the wave generation and properties. We compare Case1, Case4 and Case5 to analyze the effect of stratification on the wave generation and properties.

4.1 Generation process

Figure 2 shows distributions of horizontal velocity profiles overlaid by isopycnal surfaces at t=5, 11, 12, and 15 s in the Case1 experiment, respectively. It is seen from Fig. 2a that the turbulent wake expands uniformly in the vertical from the initial time to t=5 s, with a boundary shape of symmetric wedge type. Meanwhile, density within the wake is strongly mixing. Afterward, the wake continues expanding because the turbulent diffusive force is still stronger than the buoyancy force, and expands to its maximum extent at t=11 s (Fig. 2b). At t=12 s, wake collapse occurs and leads to the initial generation of internal waves under the inhibiting effect of buoyancy force (Fig. 2c). At t=15 s, the internal wave has escaped from the wake and begins to propagate downstream (Fig. 2d).

Figure 2 clearly illustrates the entire generation process of wake-collapse internal waves; that is, wake development, collapse and generation. The generation process basically agrees with the laboratory findings of Schooley and Stewart (1963).

4.2 Basic properties

Figure 3a-b depicts contours of density at t=20 s and 50 s in the Case1 experiment. This shows that the entire wave disturbance is almost entirely within the section between z=-110 m and z=-85 m. Zoomed views (Fig. 3c-d) of this section at t=20 s and 50 s show that the waveform in that period was essentially maintained. Moreover, displacement of isopycnal surfaces was bounded by z=-100 m, above which those surfaces were displaced upward and below which they were displaced downward. This structure indicates that the internal wave is of second-mode character (Kao and Pao, 1980; Yang et al., 2010).

To ensure that the generated internal waves are of second-mode character, consider the current field at t=50 s (Fig. 4a). The section of Fig. 4a between 375 and 390 m clearly shows symmetrical circulating cells, which coincide well with observed current (Fig. 4b) and analytical current (Fig. 4c) characteristics of a second-mode internal waves. However, given that the crucial property of internal waves is phase speed, we determined that speed numerically according to locations labeled by black dashed lines in Fig. 3, which show the wave center for comparison with the theoretical one. The theoretical phase speed was calculated via the three-layer formula of Yang et al. (2010) because the numerical result showed that the generated internal wave resembled an interfacial one in a three-layer system; the theoretical result was 0.15 m/s. The numerical phase speed was 0.17 m/s. This agreement between the theoretical and numerical phase speeds confirms that the internal wave is indeed a second-mode internal wave.

4.3 Effect of body speed 4.3.1 Effect of body speed on wave generation

A comparison among numerical results of Case1, Case2 and Case3 shows that the body speed has a strong effect on the generation of wake-collapse internal waves. In Case2, the model result did not show discernible such waves (Fig. 5). The reason, according to the laboratory results of Lin and Pao (1979), is that the buoyancy force strongly inhibits turbulent motion because of a reduction of propeller forcing in the initial stage, so that mixing in the wake is greatly reduced and there was no apparent wake collapse in the model results. Therefore, the wakecollapse internal wave was so weak that it could not be observed. This leads to the conclusion that evident wake-collapse internal waves are generated only when the speed is greater than some critical value. In Case3, the model result clearly shows that generation strength increased relative to that of Case1 and that generation time was at t=8 s (Fig. 6), 4 s earlier than in the latter case. This earlier generation is attributed to the greater speed of the moving body. Because of this faster body motion, the body is able to create more vigorous turbulent motions. Therefore, there is an earlier large-scale isopycnal overturning event near the wake boundaries, which produces an earlier wake collapse under the influence of buoyancy force.

4.3.2 Effect of body speed on wave properties

Because there was no perceptible wake-collapse internal wave in the Case2 experiment, the results from Case1 and Case3 were compared. This comparison demonstrated two notable effects of body speed on the wave properties. The first was that as body speed increased, stronger turbulent mixing caused the generation of a second-mode internal wave of larger amplitude. The wave amplitude, defined as the maximum isopycnal displacement relative to its equilibrium level, was ~1.0 m in Case1, in which the body speed was 4.5 m/s (Fig. 3b), while the wave amplitude was 1.4 m in Case3, in which the body speed was 7.5 m/s (Fig. 7b). The second effect was that the difference of body speed produced a different waveform evolution. Figure 3 shows that the waveform had almost no change from t=20 s to t=50 s in Case1, whereas it gradually evolved into a regular sinusoidal one in Case3 (Fig. 7). This suggests that the dispersion effect is stronger than the nonlinear effect in the wave generated in Case3, but this was not true of Case1. The reason is that based on internal wave theory, if such a wave's dispersion effect is much stronger than its nonlinear effect, the initial internal wave will gradually disperse toward a linear wave during its evolution. This explains the significant change of waveform in Case3. In parallel, owing to the dispersion, the horizontal length scale increased from 30 m (370-400 m in Case1) to 60 m (510-570 m in Case3). The vertical length scale was expected to decrease from 30 to 20 m from analysis of the vertical profiles of vertical velocities.

4.4 Effect of stratification magnitude 4.4.1 Effect of stratification magnitude on wave generation

A comparison among numerical results of Case1, Case4 and Case5 did not show any effect of stratification magnitude on the generation of wakecollapse internal waves. In those cases, the wake height consistently maximized at t=11 s and induced the generation of wake-collapse internal waves at t=12 s. All the generated internal waves were characterized by the second baroclinic mode.

4.4.2 Effect of stratification magnitude on wave properties

Stratification magnitude had an appreciable effect on the properties of wake-collapse internal waves. A comparison between Case1 and Case4 reveals that the internal wave generated by wake collapse in Case4 was stronger than in Case1, and also had a stronger wave-induced current. This is attributed to the weaker inhibition effect of buoyancy force in Case4. Comparison between Case4 and Case5 shows a similar result. The maximum horizontal velocity was 4.0 cm/s in Case4 (Fig. 8a) and 3.65 cm/s in Case5 (Fig. 8b). The wave amplitude decreased from 1.2 m (Case4) to 0.7 m (Case5) with increases of stratification magnitude. However, the horizontal and vertical length scales changed very little, both approximating 30 m. In short, weaker stratification produces stronger waves.

5 CONCLUSION

In this work, we constructed a numerical model that incorporates body geometry, propeller forcing, and oceanic stratification within the framework of OpenFOAM. The model was used to simulate wakecollapse internal waves under body motion at various velocities in waters with different stratification magnitudes. By analyzing the numerical results, we reached the following major conclusions.

(1) The generation and properties of wake-collapse internal waves are strongly dependent on body speed. It is only when that speed exceeds some critical value that the submerged moving body can generate such waves. That critical value is between 1.5 and 4.5 m/s. In general, if the Froude number is defined as F_r =U/ND (where N is buoyancy frequency at the depth of motion and D is the body diameter), the corresponding critical Froude number is between 14 and 26. With the increase of body speed, the wave amplitude increases and its time of generation advances. The generated wake-collapse internal waves were validated to be of second mode type. As the body speed increases, the wave's horizontal length scale increases and the waveform tends toward sinusoidal as the wave propagates;

(2) For the linearly temperature-stratified profiles, the numerical results demonstrate that the stratification magnitude has little effect on wave generation but obviously affects certain properties of wake-collapse internal waves. Specifically, the weaker the stratification, the stronger are such waves.

6 ACKNOWLEDGMENT

The work was carried out at the National Supercomputer Center in Tianjin, and calculations were performed on supercomputer TianHe-1(A).