**1. Introduction**

16 Will-be-set-by-IN-TECH

96 Computational Simulations and Applications

Tian, Z. & Dai, S. (2007). High-order compact exponential finite difference methods for

Tian, Z. & Ge, Y. (2007). A fourth-order compact ADI method for solving two-dimensional

Wood, W. (2006). An exact solution for Burgers' equation, *Communications in Numerical*

You, D. (2006). A high-order Padé ADI method for unsteady convection-diffusion equations,

Zhang, J.; Sun, H. & Zhao, J. (2002). High order compact scheme with multigrid local

2007, pp. 952 - 974.

*Mathematics*, vol. 198, 2007, pp. 268-286.

*Methods in Engineering*, Vol. 22, 2006, pp. 797 - 798.

*Journal of Computational Physics*, Vol. 214, No. 1, 2006, pp. 1-11.

*in Applied Mechanics and Engineering*, Vol. 191, 2002, pp. 4661 - 4674.

convection-diffusion type problems, *Journal of Computational Physics*, Vol. 220, No.2,

unsteady convection-diffusion problems, *Journal of Computational and Applied*

mesh refinement procedure for convection diffusion problems, *Computer Methods*

The radiation of internal gravity waves by stratified turbulent shear flows is encountered in many geophysical flows. Numerous applications include jet flows (Sutherland & Peltier 1994), grid-generated turbulence (Dohan & Sutherland 2003), boundary layers (Taylor & Sarkar 2007), collapse of mixed patches (Sutherland et al. 2007), wakes behind towed and self-propelled bodies (Lin & Pao 1979), and many others. The present paper deals with internal waves radiated by a jet flow that is created in the far wake of a sphere towed in a stratified fluid at large Froude and Reynolds numbers.

Experimental studies of internal gravity waves (IW) radiated by a towed sphere and its wake were performed both in a linearly stratified fluid by Bonneton, Chomaz & Hopfinger (1993) (further referred to as BCH) and in the presence of a thermocline by Robey (1996). The results show that internal waves radiated by the sphere (i.e. the lee waves) are stationary with respect to the sphere, and their amplitude is inversely proportional to the sphere Froude number, defined as *Fr V ND* 2 / (where *N* is the buoyancy frequency, and *V* and *D* are the towing speed and the sphere diameter). On the other hand, internal waves, radiated by the turbulent wake are not stationary with respect to the sphere, and their amplitude grows as *Fr* increases. Experimental results obtained by BCH and Robey (1996) show that if the Froude number is sufficiently large (*Fr* > 10), non-stationary IW supercede the lee waves.

Experimental results obtained by BCH show that at early times (*Nt = O*(10)) internal waves are radiated due to the collapse of the vortex coherent structures developing in the near wake (Chomaz, Bonnet & Hopfinger (1993), further referred to as CBH). The wavelength of these waves (also called "random" waves) is of the order of the sphere diameter, and their dynamics is well described theoretically under an assumption that each collapsing coherent structure can be regarded as an impulsive source of IW. Visualization of the density distribution in a horizontal plane at a distance of three sphere radii below the towing axis at times *Nt* < 40 gives a rather complicated, irregular isophase pattern of these random internal waves.

The results of BCH show also that at later times (*Nt* > 40) random internal waves are superceded by waves whose initial spatial period is of the order of five sphere diameters. The isophase distribution of these waves, although being non-stationary with respect to the sphere, is reminiscent of the regular iso-phase pattern of the lee-waves. At sufficiently late times (*Nt* > 50), the random waves disappear, and there remain only coherent IW. BCH give no explanation of the observed dynamics of these coherent internal waves.

Internal Waves Radiation by a Turbulent Jet Flow in a Stratified Fluid 99

that match parameters of a far wake of a sphere towed in a stratified fluid at large Froude and Reynolds numbers (of order *O*(10) and <sup>4</sup> *O*(10 ), respectively). The mathematical formulation and description of the numerical procedure are provided is Section 2. In Section 3 the properties of the jet flow and radiated internal waves obtained in DNS are described. In Section 4 a linear impulsive source theory is applied to explain the observed IW

A circular jet flow with the initial (reference) Gaussian profile of the streamwise (*x*) velocity

where *у* and *z* are the spanwise and vertical coordinates, respectively. Hereafter, all the variables are made dimensionless via normalizing by the velocity and length scales defined by the initial jet axis mean velocity and diameter, *U*0 and *L*<sup>0</sup> . For example, in the case of a wake created by a towed sphere, *L*0 is close to the sphere diameter ( *L D* <sup>0</sup> / 0.8 ) and *U*0 is defined by the towing speed *V* and the velocity deficit at the wake central axis at a given distance behind the sphere., e.g. *U V* <sup>0</sup> / 0.1 at *x D*/ 6 (cf. Bevilaqua & Likoudis 1978).

The Navier – Stokes equations for the fluid velocity are written under the Boussinesq

1 1 2 22 ,

Re Re

*i i ix yy ref zz ref iz*

*t i j j i ref x i ix y y ref z z ref*

*P U U U Ri*

In Eqs. (2.2) – (2.4), *Ui* ( *i xyz* , , ) and are the instantaneous deviations of the fluid velocity and density from their respective reference profiles, and *ij* is the Kroneker's symbol. The Reynolds number and the global Richardson number of the flow are defined as:

> 0 0 2 0 0

*U*

*<sup>g</sup> <sup>L</sup> Ri*

*UUUU U U U UU*

An initial linear, stable stratification of the fluid density is considered ( *ref z* ).

2 2 exp 4[ ] *U yz ref* , (2.1)

1 <sup>2</sup> RePr *t j j ref x z UU U* . (2.4)

0. *j j U* (2.3)

0 0 Re *U L* , (2.5)

, (2.6)

(2.2)

kinematics and dynamics, and final conclusions are found in Section 5.

**2. Governing equations and numerical method** 

component is considered in the form:

approximation in the dimensionless form:

The equation for the fluid density is written as:

and

The experimental study of internal waves produced by a sphere towed under a thermocline at large Froude and Reynolds numbers was performed by Robey (1996). The results also show that the turbulent wake generates coherent internal waves which are non-stationary with respect to the sphere and whose amplitude grows linearly with *Fr*. In order to explain these experimental observations Robey (1996) considers a physical model of the IW radiation by the wake. According to this model, internal waves are emitted by vortices shed by the sphere, and each vortex can be regarded as a stationary source moving at a "resonance" speed, such that its effective Froude number is of order unity. Note however, that experimental results obtained by CBH do not support the assumption that such stationary vortices are present in the near wake.

Another theoretical model of IW radiation by the wake flow has been developed by Dupont & Voisin (1996). According to this model, internal waves are emitted by an oscillating source moving with the speed of a sphere. It is assumed that such a source adequately represents coherent vortex structures which are shed by the sphere at the frequency of a spiral instability mode of the near-wake flow (observed experimentally by CBH) and collapse under gravity. Dupont & Voisin (1996) performed numerical calculations based on the Green's function formalism developed by Voisin (1994). However, numerical results of Dupont & Voisin (1996) show that the IW phase pattern predicted by the model is rather complicated and does not reproduce the regular phase distribution related to the coherent internal waves observed experimentally by BCH.

Recently there have been successful attempts to perform DNS and LES of stratified wakes by Gourlay et al. (2001), Dommermuth et al. (2002) and Brucker & Sarkar (2010). Since the mean streamwise fluid velocity in the wake is much smaller than the sphere towing speed, the variation of the flow statistical properties along the streamwise (*x*) axis can be neglected in DNS in the considered flow region. Thus the flow in DNS/LES is initialized as a sum of a circular, mean streamwise velocity profile and an *x*-periodic, random velocity component accounting for the turbulent fluctuations. This turbulent component is prescribed as a sum of Fourier harmonics with independent random phases and given amplitude power spectrum. Numerical results of Gourlay et al. (2001) (such as the temporal development of the wake axis mean velocity, width and height and the instantaneous vorticity distribution) are in good qualitative agreement with the experimental data obtained by Spedding et al. (1996) and Spedding (2001). Dommermuth et al. (2002) performed the LES study of the stratified wake and used a relaxation procedure to bring to equilibrium the initially unrelated turbulent production and dissipation. Their results are similar to the results of Gourlay et al. (2001) and show that the relaxation procedure is essential for accurately simulating the near wake but is not important if only far wake is of interest. DNS of the turbulent wake flow at high Reynolds number (Re=50000) was performed recently by Brucker & Sarkar (2010) using the same flow initialization procedure. Both Gourlay et al. (2001) and Dommermuth et al. (2002) and Brucker & Sarkar (2010) conclude that startified wakes are capable of radiating internal waves, but do not examine in detail the physical mechanism of the IW radiation and the properties of waves kinematics and dynamics. The radiation of internal waves by a turbulent jet flow in a linearly stratified fluid was studied recently by performing DNS by Druzhinin (2009). The results show that IW radiation occurs at times 10 30 *Nt* and can be described as a result of an impulsive collapse of the vertical velocity fluctuation of the initially 3D turbulent flow.

The objective of the present paper is to study in more details the process of internal waves radiation by a temporally developing turbulent jet flow in a stratified fluid by direct numerical simulation. An initially circular, turbulent jet flow is considered with parameters that match parameters of a far wake of a sphere towed in a stratified fluid at large Froude and Reynolds numbers (of order *O*(10) and <sup>4</sup> *O*(10 ), respectively). The mathematical formulation and description of the numerical procedure are provided is Section 2. In Section 3 the properties of the jet flow and radiated internal waves obtained in DNS are described. In Section 4 a linear impulsive source theory is applied to explain the observed IW kinematics and dynamics, and final conclusions are found in Section 5.

#### **2. Governing equations and numerical method**

A circular jet flow with the initial (reference) Gaussian profile of the streamwise (*x*) velocity component is considered in the form:

$$\mathcal{U}\_{ref} = \exp\left(-4\left[y^2 + z^2\right]\right),\tag{2.1}$$

where *у* and *z* are the spanwise and vertical coordinates, respectively. Hereafter, all the variables are made dimensionless via normalizing by the velocity and length scales defined by the initial jet axis mean velocity and diameter, *U*0 and *L*<sup>0</sup> . For example, in the case of a wake created by a towed sphere, *L*0 is close to the sphere diameter ( *L D* <sup>0</sup> / 0.8 ) and *U*0 is defined by the towing speed *V* and the velocity deficit at the wake central axis at a given distance behind the sphere., e.g. *U V* <sup>0</sup> / 0.1 at *x D*/ 6 (cf. Bevilaqua & Likoudis 1978). An initial linear, stable stratification of the fluid density is considered ( *ref z* ).

The Navier – Stokes equations for the fluid velocity are written under the Boussinesq approximation in the dimensionless form:

$$\begin{aligned} &\left(\hat{\sigma}\_i \mathbf{U}\_i + \mathbf{U}\_f \hat{\sigma}\_f \mathbf{U}\_i + \mathbf{U}\_{ref} \hat{\sigma}\_x \mathbf{U}\_i + \delta\_{ix} \left(\mathbf{U}\_y \hat{\sigma}\_y \mathbf{U}\_{ref} + \mathbf{U}\_z \hat{\sigma}\_z \mathbf{U}\_{ref}\right) = \\ &- \hat{\sigma}\_i P + \frac{1}{\text{Re}} \hat{\sigma}^2 \mathbf{U}\_i + \frac{1}{\text{Re}} \delta\_{ix} \left(\hat{\sigma}\_{yy}^2 \mathbf{U}\_{ref} + \hat{\sigma}\_{zz}^2 \mathbf{U}\_{ref}\right) - \delta\_{iz} \text{Ri } \mathbf{\rho}\_i \end{aligned} \tag{2.2}$$

$$
\partial\_j \mathcal{U}\_j = 0.\tag{2.3}
$$

The equation for the fluid density is written as:

$$
\hat{\sigma}\_t \mathfrak{p} + \mathcal{U}\_j \hat{\sigma}\_j \mathfrak{p} + \mathcal{U}\_{ref} \hat{\sigma}\_x \mathfrak{p} - \mathcal{U}\_z = \frac{1}{\text{Re} \, \text{Pr}} \hat{\sigma}^2 \mathfrak{p} \, . \tag{2.4}
$$

In Eqs. (2.2) – (2.4), *Ui* ( *i xyz* , , ) and are the instantaneous deviations of the fluid velocity and density from their respective reference profiles, and *ij* is the Kroneker's symbol. The Reynolds number and the global Richardson number of the flow are defined as:

$$\text{Re} = \frac{\mathcal{U}\_0 \mathcal{U}\_0}{\mathbf{v}} \,' \tag{2.5}$$

and

98 Computational Simulations and Applications

The experimental study of internal waves produced by a sphere towed under a thermocline at large Froude and Reynolds numbers was performed by Robey (1996). The results also show that the turbulent wake generates coherent internal waves which are non-stationary with respect to the sphere and whose amplitude grows linearly with *Fr*. In order to explain these experimental observations Robey (1996) considers a physical model of the IW radiation by the wake. According to this model, internal waves are emitted by vortices shed by the sphere, and each vortex can be regarded as a stationary source moving at a "resonance" speed, such that its effective Froude number is of order unity. Note however, that experimental results obtained by CBH do not support the assumption that such

Another theoretical model of IW radiation by the wake flow has been developed by Dupont & Voisin (1996). According to this model, internal waves are emitted by an oscillating source moving with the speed of a sphere. It is assumed that such a source adequately represents coherent vortex structures which are shed by the sphere at the frequency of a spiral instability mode of the near-wake flow (observed experimentally by CBH) and collapse under gravity. Dupont & Voisin (1996) performed numerical calculations based on the Green's function formalism developed by Voisin (1994). However, numerical results of Dupont & Voisin (1996) show that the IW phase pattern predicted by the model is rather complicated and does not reproduce the regular phase distribution related to the coherent

Recently there have been successful attempts to perform DNS and LES of stratified wakes by Gourlay et al. (2001), Dommermuth et al. (2002) and Brucker & Sarkar (2010). Since the mean streamwise fluid velocity in the wake is much smaller than the sphere towing speed, the variation of the flow statistical properties along the streamwise (*x*) axis can be neglected in DNS in the considered flow region. Thus the flow in DNS/LES is initialized as a sum of a circular, mean streamwise velocity profile and an *x*-periodic, random velocity component accounting for the turbulent fluctuations. This turbulent component is prescribed as a sum of Fourier harmonics with independent random phases and given amplitude power spectrum. Numerical results of Gourlay et al. (2001) (such as the temporal development of the wake axis mean velocity, width and height and the instantaneous vorticity distribution) are in good qualitative agreement with the experimental data obtained by Spedding et al. (1996) and Spedding (2001). Dommermuth et al. (2002) performed the LES study of the stratified wake and used a relaxation procedure to bring to equilibrium the initially unrelated turbulent production and dissipation. Their results are similar to the results of Gourlay et al. (2001) and show that the relaxation procedure is essential for accurately simulating the near wake but is not important if only far wake is of interest. DNS of the turbulent wake flow at high Reynolds number (Re=50000) was performed recently by Brucker & Sarkar (2010) using the same flow initialization procedure. Both Gourlay et al. (2001) and Dommermuth et al. (2002) and Brucker & Sarkar (2010) conclude that startified wakes are capable of radiating internal waves, but do not examine in detail the physical mechanism of the IW radiation and the properties of waves kinematics and dynamics. The radiation of internal waves by a turbulent jet flow in a linearly stratified fluid was studied recently by performing DNS by Druzhinin (2009). The results show that IW radiation occurs at times 10 30 *Nt* and can be described as a result of an impulsive collapse of the vertical

The objective of the present paper is to study in more details the process of internal waves radiation by a temporally developing turbulent jet flow in a stratified fluid by direct numerical simulation. An initially circular, turbulent jet flow is considered with parameters

stationary vortices are present in the near wake.

internal waves observed experimentally by BCH.

velocity fluctuation of the initially 3D turbulent flow.

$$Ri = \frac{g\Delta\rho\_0}{\rho\_0} \frac{L\_0}{\mathcal{U}\_0^2} \, ^\prime \tag{2.6}$$

Internal Waves Radiation by a Turbulent Jet Flow in a Stratified Fluid 101

development and IW radiation in a two-dimensional jet flow. In this approach, the absorption of internal waves radiated from the jet core is ensured by increasing the fluid kinematic viscosity in a layer near the boundaries of the computational domain. A similar viscous-sponge boundary condition was successfully employed recently by Taylor & Sarkar (2007) in their LES study of internal waves radiated by a turbulent boundary layer. Thus, in our DNS procedure the viscosity is increased by a factor of 20 in the region 1/2 2 2 *y z* <sup>7</sup> ,

as compared to the viscosity in the rest of the computational domain (where it is constant and equals 1/Re). Our DNS results show that such artificial viscosity enhancement does not affect the jet flow dynamics and ensures complete absorption of internal waves radiated

In order to define the flow statistical characteristics at a given time moment (*t*), the *x* –

*i i ix ref*

*U y z t U x y z t dx U y z <sup>L</sup>*

<sup>1</sup> ( , ,) ( , , ,) ( , )

0 <sup>1</sup> ( , ,) ( , , ,) *Lx*

*x*

where *i xyz* , , , and 36 *Lx* . The jet maximum mean velocity (*Um* ) and its width and

<sup>1</sup> ( ) ( , 0, ) *y x*

<sup>1</sup> ( ) ( 0, , ) *z x*

*L t U y z t dz <sup>U</sup>*

The initial, random component of the fluid velocity is prescribed as a 3D homogeneous isotropic turbulence field consisting of a sum of Fourier harmonics with independent random phases and an isotropic amplitude power spectrum (Gourlay et al. 2001). We prescribe the initial turbulence spectrum to be uniform over the wavenumber range 0.5 5 *k* and equal to zero for other *k*'s (Fig. 1). Thus, the initial spectrum includes the

*L t U y z t dy <sup>U</sup>*

*m*

*m*

The instantaneous r.m.s. fluctuations of the velocity and density fields,

and their amplitudes (the respective absolute maxima) are also evaluated.

, (2.10)

*<sup>y</sup> z t x y z t dx <sup>L</sup>* , (2.11)

*Ut U m x* ( ) max ( , , ) *y z t* , (2.12)

, (2.13)

. (2.14)

2 21/2 '( , ,) *U yzt U U i ii* (2.15)

2 2 1/2 '( , , ) *yzt* , (2.16)

0

*x*

*Lx*

from the jet core region toward the boundaries.

height ( *Ly*,*<sup>z</sup>* ) are further evaluated as:

and

averaged velocity and density fields are evaluated as:

where is the fluid kinematic viscosity, *g* the acceleration due to gravity, and ( 0 0 / ) the absolute relative density variation in the vertical direction on scale *L*<sup>0</sup> . (Note that since *U V* <sup>0</sup> / 0.1 , the Reynolds and Froude numbers of the sphere are about 10 times larger than Re and 1/2 *Fr Ri* , respectively.) As follows from eqs. (2.2) - (2.6), the dimensionless buoyancy frequency equals 1/2 *N Ri Fr* 1 / . The Prandtl number is set equal to unity. The results of DNS by Stadler et al. (2010) performed recently show that that Pr = 1 is a reasonable approximation for Pr = 7 which is the case of the thermal stratification in the ocean.

In order to perform numerical simulation in a region finite over the vertical coordinate and provide a sufficient resolution of the flow field, a mapping for the *z*-coordinate is employed in the form:

tanh 9 *<sup>z</sup>* , (2.7)

so that

$$z = 4.5 \ln \left( \frac{1 + \xi}{1 - \xi} \right) \tag{2.8}$$

and 1 1 for *z* . The derivatives over *z*-coordinate in Eqs. (2.2)-(2.4) are rewritten as:

$$\frac{\partial}{\partial z} = \left(\frac{1-\xi^2}{9}\right)\frac{\partial}{\partial \xi}\,. \tag{2.9}$$

Equations (2.2) - (2.4) are discretized in a parallelepiped domain with sizes 0 36 *x* , 9 9 *y* and 1 1 by employing a finite difference method of the secondorder accuracy on a uniform rectangular staggered grid consisting of 480 × 240 × 240 nodes in *x, y*, and *ξ*(*z*) directions, respectively. The terms on the left hand side of equations (2.2) and (2.4), *U U ref x i* and *Uref <sup>x</sup>* , that are responsible for the advection by the reference velocity ( , ) *Uref y z* , are evaluated at each time step by Fourier interpolation which is a more accurate method to compute these terms as compared to a pure finite difference scheme (Gerz et al. (1989)). The integration is advanced in time using the Adams-Bashforth method with time step *t* = 0.0075. The shear-free (Neumann) boundary condition is prescribed in the spanwise and vertical directions (at y = ± 9 and = ± 1), and *x*-periodic boundary condition is prescribed in the streamwise direction. The Poisson equation for the pressure is solved by fast Fourier transform over x-coordinate, cosine transform over y-coordinate, and Gauss elimination method over z-coordinate.

As it was observed by Gourlay et al. (2001) in their DNS study, internal waves radiated by the jet turbulence propagate away from the jet core toward the boundaries of the computational domain in the *y* and *z* directions. Thus one needs to avoid reflection of the waves from the boundaries and their coming back into the flow region. One way to do this is to prescribe the Sommerfeld radiation boundary condition (cf. e.g. Javam et al. 2000). In the present study we adopt a less computationally expensive approach analogous to that employed by Sutherland & Peltier (1994) in their numerical study of the instability development and IW radiation in a two-dimensional jet flow. In this approach, the absorption of internal waves radiated from the jet core is ensured by increasing the fluid kinematic viscosity in a layer near the boundaries of the computational domain. A similar viscous-sponge boundary condition was successfully employed recently by Taylor & Sarkar (2007) in their LES study of internal waves radiated by a turbulent boundary layer. Thus, in

our DNS procedure the viscosity is increased by a factor of 20 in the region 1/2 2 2 *y z* <sup>7</sup> ,

as compared to the viscosity in the rest of the computational domain (where it is constant and equals 1/Re). Our DNS results show that such artificial viscosity enhancement does not affect the jet flow dynamics and ensures complete absorption of internal waves radiated from the jet core region toward the boundaries.

In order to define the flow statistical characteristics at a given time moment (*t*), the *x* – averaged velocity and density fields are evaluated as:

$$<\mathcal{U}\_i > \{y, z, t\} = \frac{1}{L\_\chi} \int\_0^L \mathcal{U}\_i(x, y, z, t) dx + \mathcal{S}\_{ix} \mathcal{U}\_{ref}(y, z) \,. \tag{2.10}$$

$$<\rho>(y,z,t) = \frac{1}{L\_x} \int\_0^{L\_y} \rho(x,y,z,t)dx\,\,\,\tag{2.11}$$

where *i xyz* , , , and 36 *Lx* . The jet maximum mean velocity (*Um* ) and its width and height ( *Ly*,*<sup>z</sup>* ) are further evaluated as:

$$\mathcal{U}I\_m(t) = \max\left\{<\mathcal{U}\_x>(y,z,t)\right\},\tag{2.12}$$

$$L\_y(t) = \frac{1}{\mathcal{U}\_m} \int\_{-\infty}^{\infty} \zeta \mathcal{U}\_x > (y, z = 0, t) dy \,\prime,\tag{2.13}$$

$$L\_z(t) = \frac{1}{\mathcal{U}\_m} \int\_{-\infty}^{\infty} \zeta \mathcal{U}\_x > (y = 0, z, t) dz \,. \tag{2.14}$$

The instantaneous r.m.s. fluctuations of the velocity and density fields,

$$\mathcal{U}I\_i'(y, z, t) = <\mathcal{U}\_i^2 - <\mathcal{U}\_i >^2 ^{1/2} \tag{2.15}$$

and

100 Computational Simulations and Applications

where is the fluid kinematic viscosity, *g* the acceleration due to gravity, and ( 0 0 / ) the absolute relative density variation in the vertical direction on scale *L*<sup>0</sup> . (Note that since *U V* <sup>0</sup> / 0.1 , the Reynolds and Froude numbers of the sphere are about 10 times larger than Re and 1/2 *Fr Ri* , respectively.) As follows from eqs. (2.2) - (2.6), the dimensionless buoyancy frequency equals 1/2 *N Ri Fr* 1 / . The Prandtl number is set equal to unity. The results of DNS by Stadler et al. (2010) performed recently show that that Pr = 1 is a reasonable approximation for Pr = 7 which is the case of the thermal stratification in the

In order to perform numerical simulation in a region finite over the vertical coordinate and provide a sufficient resolution of the flow field, a mapping for the *z*-coordinate is employed

> tanh 9

<sup>1</sup> 4.5ln 1

and 1 1 for *z* . The derivatives over *z*-coordinate in Eqs. (2.2)-(2.4) are

<sup>2</sup> 1

Equations (2.2) - (2.4) are discretized in a parallelepiped domain with sizes 0 36 *x* , 9 9 *y* and 1 1 by employing a finite difference method of the secondorder accuracy on a uniform rectangular staggered grid consisting of 480 × 240 × 240 nodes in *x, y*, and *ξ*(*z*) directions, respectively. The terms on the left hand side of equations (2.2) and (2.4), *U U ref x i* and *Uref <sup>x</sup>* , that are responsible for the advection by the reference velocity ( , ) *Uref y z* , are evaluated at each time step by Fourier interpolation which is a more accurate method to compute these terms as compared to a pure finite difference scheme (Gerz et al. (1989)). The integration is advanced in time using the Adams-Bashforth method with time step *t* = 0.0075. The shear-free (Neumann) boundary condition is prescribed in the spanwise and vertical directions (at y = ± 9 and = ± 1), and *x*-periodic boundary condition is prescribed in the streamwise direction. The Poisson equation for the pressure is solved by fast Fourier transform over x-coordinate, cosine transform over y-coordinate, and

As it was observed by Gourlay et al. (2001) in their DNS study, internal waves radiated by the jet turbulence propagate away from the jet core toward the boundaries of the computational domain in the *y* and *z* directions. Thus one needs to avoid reflection of the waves from the boundaries and their coming back into the flow region. One way to do this is to prescribe the Sommerfeld radiation boundary condition (cf. e.g. Javam et al. 2000). In the present study we adopt a less computationally expensive approach analogous to that employed by Sutherland & Peltier (1994) in their numerical study of the instability

*z* 9 

*z*

*<sup>z</sup>* , (2.7)

, (2.8)

. (2.9)

ocean.

so that

rewritten as:

Gauss elimination method over z-coordinate.

in the form:

$$\mathfrak{p}'(y, z, t) = <\mathfrak{p}^2 - <\mathfrak{p} >^2 >^{1/2},\tag{2.16}$$

and their amplitudes (the respective absolute maxima) are also evaluated. The initial, random component of the fluid velocity is prescribed as a 3D homogeneous isotropic turbulence field consisting of a sum of Fourier harmonics with independent random phases and an isotropic amplitude power spectrum (Gourlay et al. 2001). We prescribe the initial turbulence spectrum to be uniform over the wavenumber range 0.5 5 *k* and equal to zero for other *k*'s (Fig. 1). Thus, the initial spectrum includes the

Internal Waves Radiation by a Turbulent Jet Flow in a Stratified Fluid 103

U L L

y z

> t 1/3

t -2/3

L , Ri = 1 L , Ri = 1 L , Ri = 3

y

y z

z

t 0.33

0.1 1 10 100

0.1 1 10

Fig. 2. Temporal development of the velocity and length scales of the non-stratified jet flow.

In order to verify the performance of our numerical procedure we performed DNS with Re = 700 of a homogeneous jet flow (where Ri = 0) and a stratified jet flow with two different Richardson numbers (Ri = 1 and Ri = 3) and compared our numerical results with well-

0.1

Nt Nt

Fig. 3. Temporal development of the jet mean velocity maximum *Um* (left) and its width *Ly* and height *Lz* (in dashed line) (right). Asymptotics 0.75 *U t <sup>m</sup>* <sup>~</sup> and 0.33 *L t y z*, ~ are shown

1

10

L , Ri = 3 Um Ly,z

t

0.1

known asymptotics of homogeneous and stratified far wakes.

U , Ri=1 U , Ri=3

m m

0.1 1 10 100

t -0.75

**3. Numerical results 3.1 Jet flow properties** 

0.1

in dotted line.

1

1

10

most unstable (spiral) mode with *k* = 1.46 of an inviscid non-stratified circular jet flow (Batchelor and Gill 1962), as well as a higher-*k* axisymmetrical Kelvin-Helmholtz instability mode similar to that observed in the near wake flow (CBH 1993).

At present it is an open question whether the spiral shedding mode observed in the near wake (CBH 1993) is related to the spiral instability mode of the *x*-periodic circular jet, although the Strouhal numbers of the two modes are close (for the spiral mode St = *k*/2 = 0.2 whereas for the shedding spiral mode St = 0.175, CBH 1993). It is possible that the two modes are related due to some correlation between the development of the instability of the near-wake flow and the location of the separation point of the boundary layer at the sphere and the development of the spiral shedding mode. It is important to note however that the development of the spiral instability mode is an intrinsic property of the circular jet/wake flow itself and not related to presence of the sphere.

Since at present there is no knowledge concerning the shape of the spectrum in the homogeneous near wake, we prescribe an initial uniform-amplitude spectrum. The initial amplitude of the random velocity component is prescribed to be about 40% of the mean velocity maximum. This value is close to the turbulence intensity measured by Bevilaqua and Lykoudis (1978) in a homogeneous turbulent wake at the distance / 6 *x D* behind the sphere. The initial velocity field is further windowed with the reference profile *zyU* ),( *ref*

(2.1) and made divergence-free , similarly to the initialization procedure employed by Gourlay et al. (2001).

The initial deviation of the density from the reference linear profile is prescribed to be zero. This is in agreement with the experimental observations of BCH showing that the effect of the turbulent mixing in the near wake of a towed sphere is weak and its contribution to the radiated IW field is negligible.

Fig. 1. The initial velocity spectrum.

Fig. 2. Temporal development of the velocity and length scales of the non-stratified jet flow.

### **3. Numerical results**

102 Computational Simulations and Applications

most unstable (spiral) mode with *k* = 1.46 of an inviscid non-stratified circular jet flow (Batchelor and Gill 1962), as well as a higher-*k* axisymmetrical Kelvin-Helmholtz instability

At present it is an open question whether the spiral shedding mode observed in the near wake (CBH 1993) is related to the spiral instability mode of the *x*-periodic circular jet, although the Strouhal numbers of the two modes are close (for the spiral mode St = *k*/2 = 0.2 whereas for the shedding spiral mode St = 0.175, CBH 1993). It is possible that the two modes are related due to some correlation between the development of the instability of the near-wake flow and the location of the separation point of the boundary layer at the sphere and the development of the spiral shedding mode. It is important to note however that the development of the spiral instability mode is an intrinsic property of the circular jet/wake

Since at present there is no knowledge concerning the shape of the spectrum in the homogeneous near wake, we prescribe an initial uniform-amplitude spectrum. The initial amplitude of the random velocity component is prescribed to be about 40% of the mean velocity maximum. This value is close to the turbulence intensity measured by Bevilaqua and Lykoudis (1978) in a homogeneous turbulent wake at the distance / 6 *x D* behind the sphere. The initial velocity field is further windowed with the reference profile *zyU* ),( *ref* (2.1) and made divergence-free , similarly to the initialization procedure employed by

The initial deviation of the density from the reference linear profile is prescribed to be zero. This is in agreement with the experimental observations of BCH showing that the effect of the turbulent mixing in the near wake of a towed sphere is weak and its contribution to the

> U U U

y x z

0.1 1 10

k

mode similar to that observed in the near wake flow (CBH 1993).

flow itself and not related to presence of the sphere.

0.001

Fig. 1. The initial velocity spectrum.

0.01

0.1

1

Gourlay et al. (2001).

radiated IW field is negligible.

#### **3.1 Jet flow properties**

In order to verify the performance of our numerical procedure we performed DNS with Re = 700 of a homogeneous jet flow (where Ri = 0) and a stratified jet flow with two different Richardson numbers (Ri = 1 and Ri = 3) and compared our numerical results with wellknown asymptotics of homogeneous and stratified far wakes.

Fig. 3. Temporal development of the jet mean velocity maximum *Um* (left) and its width *Ly* and height *Lz* (in dashed line) (right). Asymptotics 0.75 *U t <sup>m</sup>* <sup>~</sup> and 0.33 *L t y z*, ~ are shown in dotted line.

Internal Waves Radiation by a Turbulent Jet Flow in a Stratified Fluid 105

The density fluctuations amplitude ' increases and reaches its maximum at 2 *Nt* due to the transfer of the initial turbulence kinetic energy into the potential energy of the displaced isopycnals. Then ' decreases, and this reduction is accompanied by the growth of the vertical velocity fluctuations until 3.5 *Nt* . At later times the anti-phase oscillations of ' and ' *U <sup>z</sup>* are smeared out due to their decay. Similar behavior was observed by Dommermuth et al. (2002) and interpreted as the radiation of internal waves by collapsing

0.1 1 10 100

k

0.1 1 10 100

Nt=90

1e-007 1e-006 1e-005 0.0001 0.001 0.01 0.1

Fig. 5. Velocity fluctuations spatial spectra at different time moments. Asymptotics 3/5 *k* is

Figure 5 shows the velocity wavenumber spectra ( ) *U k <sup>i</sup>* of the *i* = *x*, *y*, *z* velocity components obtained in DNS for Ri = 1 at different time moments. Each spectrum is evaluated by a

1e-005

Nt=30

0.0001

0.001

0.01

0.1

U

x

y

z

U

U

Nt=3 Nt=6

0.1 1 10 100

k

shown in dotted line. Here and below in Figs. 6, 7 Ri = 1.

0.1 1 10 100

1e-006

1e-005

0.0001

0.001

0.01

0.1

1e-005

0.0001

0.001

0.01

0.1

vertical velocity fluctuations in the wake.

Figures 2 and 3 compare the temporal development of the flow integral scales, such as the mean velocity maximum, and jet width and height obtained in our DNS, with the asymptotics 2/3 *U t <sup>m</sup>* <sup>~</sup> and 0.75 *U t <sup>m</sup>* <sup>~</sup> and 0.33 *L t y z*, <sup>~</sup> , observed in homogeneous and stratified wakes in laboratory experiments (Bevilaqua & Lykoudis 1978, Spedding 1997). In the case of a homogeneous wake, the agreement is very good for the considered times *t* 20 (cf. Fig. 2). (It is important to note that since we normalize by the velocity deficit *U V* <sup>0</sup> / 0.1 , time moment *t* = 20 corresponds to the distance / / 200 <sup>0</sup> *x D tV U* behind the sphere.) In the stratified case, the agreement with the asymptotics 0.75 *U t <sup>m</sup>* <sup>~</sup> and 0.33 *L t y z*, ~ is also quite satisfactory.

Figure 4 presents the temporal development of the velocity and density fluctuations amplitudes for Ri = 1 and 3. The figure shows that at early times ( 10 *Nt* ) the velocity fluctuations amplitudes are of the same order, i.e. the flow remains nearly isotropic. At later times, the vertical velocity decreases much faster as compared to the horizontal velocity components and at times 30 *Nt* the amplitude ' *U <sup>z</sup>* becomes almost by two orders of magnitude smaller than the horizontal velocity amplitudes ' *U <sup>x</sup>* and ' *U <sup>y</sup>* . Thus, at that stage, the flow becomes quasi-two-dimensional.

Fig. 4. Temporal development of the velocity and density fluctuations amplitudes for Ri = 1 (left) and Ri = 3 (right).

Figures 2 and 3 compare the temporal development of the flow integral scales, such as the mean velocity maximum, and jet width and height obtained in our DNS, with the asymptotics 2/3 *U t <sup>m</sup>* <sup>~</sup> and 0.75 *U t <sup>m</sup>* <sup>~</sup> and 0.33 *L t y z*, <sup>~</sup> , observed in homogeneous and stratified wakes in laboratory experiments (Bevilaqua & Lykoudis 1978, Spedding 1997). In the case of a homogeneous wake, the agreement is very good for the considered times *t* 20 (cf. Fig. 2). (It is important to note that since we normalize by the velocity deficit *U V* <sup>0</sup> / 0.1 , time moment *t* = 20 corresponds to the distance / / 200 <sup>0</sup> *x D tV U* behind the sphere.) In the stratified case, the agreement with the asymptotics 0.75 *U t <sup>m</sup>* <sup>~</sup> and

Figure 4 presents the temporal development of the velocity and density fluctuations amplitudes for Ri = 1 and 3. The figure shows that at early times ( 10 *Nt* ) the velocity fluctuations amplitudes are of the same order, i.e. the flow remains nearly isotropic. At later times, the vertical velocity decreases much faster as compared to the horizontal velocity components and at times 30 *Nt* the amplitude ' *U <sup>z</sup>* becomes almost by two orders of magnitude smaller than the horizontal velocity amplitudes ' *U <sup>x</sup>* and ' *U <sup>y</sup>* . Thus, at that

> Ux' Uy' Uz'

0.0001

Fig. 4. Temporal development of the velocity and density fluctuations amplitudes for Ri = 1

0.001

0.01

0.1

1

Ri = 1 Ri = 3

0.1 1 10 100

Nt

0.33 *L t y z*, ~ is also quite satisfactory.

stage, the flow becomes quasi-two-dimensional.

0.1 1 10 100

Nt

0.0001

(left) and Ri = 3 (right).

0.001

0.01

0.1

1

The density fluctuations amplitude ' increases and reaches its maximum at 2 *Nt* due to the transfer of the initial turbulence kinetic energy into the potential energy of the displaced isopycnals. Then ' decreases, and this reduction is accompanied by the growth of the vertical velocity fluctuations until 3.5 *Nt* . At later times the anti-phase oscillations of ' and ' *U <sup>z</sup>* are smeared out due to their decay. Similar behavior was observed by Dommermuth et al. (2002) and interpreted as the radiation of internal waves by collapsing vertical velocity fluctuations in the wake.

Fig. 5. Velocity fluctuations spatial spectra at different time moments. Asymptotics 3/5 *k* is shown in dotted line. Here and below in Figs. 6, 7 Ri = 1.

Figure 5 shows the velocity wavenumber spectra ( ) *U k <sup>i</sup>* of the *i* = *x*, *y*, *z* velocity components obtained in DNS for Ri = 1 at different time moments. Each spectrum is evaluated by a

Internal Waves Radiation by a Turbulent Jet Flow in a Stratified Fluid 107

x 0 5 10 15 20 25 30 35

central (*x,y*) plane at time moments *Nt* = 6 (top frame) and *Nt* = 90 (bottom frame). The

Figure 7 shows the instantaneous distributions of the r.m.s. of the density and horizontal

(2.15) and (2.16) at different time moments. The figure shows that at sufficiently late times (*Nt* > 10) the vertical velocity fluctuations, '( , ) *U yz <sup>z</sup>* , are mostly present in the region | | 1 *z* whereas the horizontal component, '( , ) *U yz <sup>h</sup>* , is most pronounced in a horizontal layer ||1 *z* in the region | | *<sup>y</sup> y L* . At times *Nt* = 15, 30 the vertical velocity fluctuations are most pronounced in the elongated regions oriented at the angle *θ* with respect to the vertical axis in the range 0 0 40 60 . It is a well-known property of stratified turbulence that at sufficiently late times, after the initial collapse of 3D turbulence under the action of buoyancy force, the vertical velocity is associated mostly with propagating internal waves (IW) whereas the horizontal velocity is associated with quasi-2D fluid motions. Fig. 7 shows that at later times (*Nt* = 60, 90) amplitude ' *U <sup>z</sup>* decreases by the order of magnitude as compared to the time moment *Nt* = 15, and the regions of local maxima of '( , ) *U yz <sup>z</sup>* are getting aligned with the vertical axis. This picture is similar to the one predicted by a linear theory for the IW radiated by an impulsive point source (Zavolski & Zaytsev 1984, Voisin 1991). According to this theory, amplitude of IW emitted impulsively by a point source has its maximum along the iso-phase line oriented at <sup>0</sup> arctan 2 55 with respect to the

1/2 2 2 '( , ) ' ' *U yz U U h xy* and '( , ) *U yz <sup>z</sup>* , evaluated from

*<sup>z</sup>* in the horizontal

Fig. 6. Instantaneous contours of the vertical vorticity component

and vertical velocity fields, '( , ) *y z* ,

mean flow is from left to right along the *x-*axis.

**3.2 Internal waves radiation** 


z


y


0

5


0

5

Fourier transform of the velocity *i*-component (,,) *U xyz <sup>i</sup>* over the streamwise (*x*) coordinate and averaging over 9 *x*-realizations located in the jet core region 2 2 *y z* 0.15 . The figure shows that at early time moments *Nt* = 3, 6 the spectra ( ) *U k <sup>i</sup>* are of the same magnitude, i.e. the velocity field remains isotropic. It is important to note that at these early times the spectra are characterized by distinctive peaks at *<sup>s</sup> k k* (for *y*-component) and 0.5 *<sup>s</sup> k k* (for *x*-component) where 1.46 *<sup>s</sup> k* is the wavenumber of the spiral instability mode of the nonstratified jet flow (Batchelor & Gill 1962). At later times (*Nt* =30, 90) most of the energy becomes accumulated in the horizontal velocity components ( ) *U k <sup>x</sup>* and ( ) *U k <sup>y</sup>* with pronounced peaks at 0.5 *<sup>s</sup> k k* .

The dynamics of the spectrum in Fig. 5 can be interpreted as follows. At early times (*Nt* < 10) there occurs a preferential growth of the most unstable mode which, in a certain sense, is analogous to the spiral instability mode of the non-stratified jet flow (Batchelor & Gill 1962). (This is in agreement with numerous experimental observations (cf. e.g. Spedding et al. (1996)) showing that at early times the stratified wake flow develops similarly to the homogeneous wake flow.) It is important to note that this instability mode is an *intrinsic* property of the jet flow itself and is not directly related to the shedding instability mode of the near-wake flow observed experimentally by CBH (1993). It is important to note also that, in our DNS, no selected modes are initially present (cf. Fig. 1). At later times (*Nt* > 30) the flow is governed by the non-linear interaction and competition between the most developed flow modes which leads to the preferential growth of a quasi-2D mode with wavenumber 0.5 *<sup>s</sup> k* .

Figure 5 shows also that the fluid kinetic energy is negligible for wavenumbers larger than max *k* 20 (the corresponding spatial scale is min max *l k* 2 / 0.3 ). Thus, for the considered spatial resolution (Δ*x* = 0.075) the flow, including the smallest scales, is well-resolved in DNS.

Figure 6 presents the instantaneous distribution of the flow vorticity *z-*component, ( ) *z xU UU y y x ref* in the plane (*x, y*) at *z* = 0 obtained in DNS for Ri = 1 at time moments *Nt* = 6 and *Nt* = 90. The figure shows that at late times the vorticity distribution is characterized by the development of large-scale vortices of alternating polarity arranged in the horizontal plane (*x, y*) in the vicinity of the jet streamwise axis. Similar vortex structures (also called the "pancake" eddies) are known to be the common feature of stratified far wakes and observed both in laboratory experiments (cf. e.g. Spedding et al. (1996)) and in numerical simulations (Gourlay et al. (2001), Dommermuth et al. (2002)).

It is of interest to note that at times 30 < *Nt* < 90 the amplitude ' *U <sup>y</sup>* of the spanwise fluctuation velocity component remains almost constant whereas the streamwise velocity amplitude, ' *U <sup>x</sup>* , decreases with time, so that at sufficiently late times (*Nt* > 60) the spanwise velocity fluctuations prevail over the streamwise fluctuations (cf. Fig. 4). This behavior of ' *U <sup>y</sup>* and ' *U <sup>x</sup>* is more pronounced for a larger Richardson number (cf. cases Ri= 1 and Ri = 3 in Fig. 4). Figure 5 shows also that at that stage a distinct peak develops in the horizontal velocity spectra , ( ) *U k x y* at wavenumber 0.5 *<sup>s</sup> k k* . This behavior of the spanwise velocity amplitude and the spectrum is in general agreement with the hypothesis put forward by Spedding (2001) who pointed out that an *intrinsic* mean flow sinuous instability may be responsible for the development of the large-scale coherent vortex structures in late stratified wakes.

Fourier transform of the velocity *i*-component (,,) *U xyz <sup>i</sup>* over the streamwise (*x*) coordinate and averaging over 9 *x*-realizations located in the jet core region 2 2 *y z* 0.15 . The figure shows that at early time moments *Nt* = 3, 6 the spectra ( ) *U k <sup>i</sup>* are of the same magnitude, i.e. the velocity field remains isotropic. It is important to note that at these early times the spectra are characterized by distinctive peaks at *<sup>s</sup> k k* (for *y*-component) and 0.5 *<sup>s</sup> k k* (for *x*-component) where 1.46 *<sup>s</sup> k* is the wavenumber of the spiral instability mode of the nonstratified jet flow (Batchelor & Gill 1962). At later times (*Nt* =30, 90) most of the energy becomes accumulated in the horizontal velocity components ( ) *U k <sup>x</sup>* and ( ) *U k <sup>y</sup>* with

The dynamics of the spectrum in Fig. 5 can be interpreted as follows. At early times (*Nt* < 10) there occurs a preferential growth of the most unstable mode which, in a certain sense, is analogous to the spiral instability mode of the non-stratified jet flow (Batchelor & Gill 1962). (This is in agreement with numerous experimental observations (cf. e.g. Spedding et al. (1996)) showing that at early times the stratified wake flow develops similarly to the homogeneous wake flow.) It is important to note that this instability mode is an *intrinsic* property of the jet flow itself and is not directly related to the shedding instability mode of the near-wake flow observed experimentally by CBH (1993). It is important to note also that, in our DNS, no selected modes are initially present (cf. Fig. 1). At later times (*Nt* > 30) the flow is governed by the non-linear interaction and competition between the most developed flow modes which leads to the preferential growth of a quasi-2D mode with

Figure 5 shows also that the fluid kinetic energy is negligible for wavenumbers larger than max *k* 20 (the corresponding spatial scale is min max *l k* 2 / 0.3 ). Thus, for the considered spatial resolution (Δ*x* = 0.075) the flow, including the smallest scales, is well-resolved in

Figure 6 presents the instantaneous distribution of the flow vorticity *z-*component, ( ) *z xU UU y y x ref* in the plane (*x, y*) at *z* = 0 obtained in DNS for Ri = 1 at time moments *Nt* = 6 and *Nt* = 90. The figure shows that at late times the vorticity distribution is characterized by the development of large-scale vortices of alternating polarity arranged in the horizontal plane (*x, y*) in the vicinity of the jet streamwise axis. Similar vortex structures (also called the "pancake" eddies) are known to be the common feature of stratified far wakes and observed both in laboratory experiments (cf. e.g. Spedding et al. (1996)) and in

It is of interest to note that at times 30 < *Nt* < 90 the amplitude ' *U <sup>y</sup>* of the spanwise fluctuation velocity component remains almost constant whereas the streamwise velocity amplitude, ' *U <sup>x</sup>* , decreases with time, so that at sufficiently late times (*Nt* > 60) the spanwise velocity fluctuations prevail over the streamwise fluctuations (cf. Fig. 4). This behavior of

 *U <sup>y</sup>* and ' *U <sup>x</sup>* is more pronounced for a larger Richardson number (cf. cases Ri= 1 and Ri = 3 in Fig. 4). Figure 5 shows also that at that stage a distinct peak develops in the horizontal velocity spectra , ( ) *U k x y* at wavenumber 0.5 *<sup>s</sup> k k* . This behavior of the spanwise velocity amplitude and the spectrum is in general agreement with the hypothesis put forward by Spedding (2001) who pointed out that an *intrinsic* mean flow sinuous instability may be responsible for the development of the large-scale coherent vortex structures in late

numerical simulations (Gourlay et al. (2001), Dommermuth et al. (2002)).

pronounced peaks at 0.5 *<sup>s</sup> k k* .

wavenumber 0.5 *<sup>s</sup> k* .

DNS.

'

stratified wakes.

Fig. 6. Instantaneous contours of the vertical vorticity component *<sup>z</sup>* in the horizontal central (*x,y*) plane at time moments *Nt* = 6 (top frame) and *Nt* = 90 (bottom frame). The mean flow is from left to right along the *x-*axis.

#### **3.2 Internal waves radiation**

Figure 7 shows the instantaneous distributions of the r.m.s. of the density and horizontal and vertical velocity fields, '( , ) *y z* , 1/2 2 2 '( , ) ' ' *U yz U U h xy* and '( , ) *U yz <sup>z</sup>* , evaluated from (2.15) and (2.16) at different time moments. The figure shows that at sufficiently late times (*Nt* > 10) the vertical velocity fluctuations, '( , ) *U yz <sup>z</sup>* , are mostly present in the region | | 1 *z* whereas the horizontal component, '( , ) *U yz <sup>h</sup>* , is most pronounced in a horizontal layer ||1 *z* in the region | | *<sup>y</sup> y L* . At times *Nt* = 15, 30 the vertical velocity fluctuations are most pronounced in the elongated regions oriented at the angle *θ* with respect to the vertical axis in the range 0 0 40 60 . It is a well-known property of stratified turbulence that at sufficiently late times, after the initial collapse of 3D turbulence under the action of buoyancy force, the vertical velocity is associated mostly with propagating internal waves (IW) whereas the horizontal velocity is associated with quasi-2D fluid motions. Fig. 7 shows that at later times (*Nt* = 60, 90) amplitude ' *U <sup>z</sup>* decreases by the order of magnitude as compared to the time moment *Nt* = 15, and the regions of local maxima of '( , ) *U yz <sup>z</sup>* are getting aligned with the vertical axis. This picture is similar to the one predicted by a linear theory for the IW radiated by an impulsive point source (Zavolski & Zaytsev 1984, Voisin 1991). According to this theory, amplitude of IW emitted impulsively by a point source has its maximum along the iso-phase line oriented at <sup>0</sup> arctan 2 55 with respect to the

Internal Waves Radiation by a Turbulent Jet Flow in a Stratified Fluid 109

Figure 8 shows the instantaneous distribution of the vertical velocity component, *Uz* , in the vertical (*x, z*) planes at *y* = 0 and *y* = 1.5 obtained in DNS at time moments *Nt* = 15 and *Nt* = 90. The figure shows that at time moment *Nt* = 15, in the region sufficiently far away from the jet streamwise axis (for *z* 1 ), IW isophase lines, which are associated with the regions of maxima and minima of *Uz* , are mostly in the direction of the mean flow and oriented at the angle of <sup>0</sup> 50 with respect to the vertical. The velocity distribution has a nearly periodical spatial structure with a period of about 5 dimensionless units along the *x*coordinate. The figure shows also that at time *Nt* = 90 the vertical velocity amplitude is diminished by the order of magnitude and its maxima and minima are arranged in columns with the same spatial periodicity. These maxima and minima can be associated with smallamplitude buoyancy oscillations which supercede impulsively radiated IW at that stage of

Uz



X

Fig. 8. Instantaneous contours of the vertical velocity *U <sup>z</sup>* in the vertical (*x,z*) plane at *y* = 0 (left column) and *y* = 1.5 (right column) at time moments *Nt* = 15 (top frame) and *Nt* = 90

Figure 9 presents the distribution of the vertical velocity, *Uz* , in the horizontal (*x, y*)-plane at the distance *z* = 3 above the jet streamwise axis obtained in DNS at time moments *Nt* = 12, 15, 18. The figure shows that the IW patterns in that plane are similar to the IW patterns of stationary lee-waves and characterized by a well-defined streamwise spatial period of about 5 dimensionless units (or about 4-5 sphere diameters in the dimensional form). A similar isophase pattern was observed by BCH in the laboratory experiment in a horizontal plane at a distance of three sphere radii away from the towing axis. Note also that this spatial period is close to the wavelength of the spiral instability mode of the non-stratified jet flow

Figure 10 shows the temporal development of the vertical velocity and density, *Uz* and *ρ*, obtained in DNS at the point located at *x* = 18, with the angle with respect to the vertical *θ* =

0 5 10 15 20 25 30 35

the flow development (Voisin 1991).

0 5 10 15 20 25 30 35

2 / 4.3 *s s k* for 1.46 *<sup>s</sup> k* (Batchelor and Gill 1962).

Z



(bottom frame).

0

5

0

5

vertical. At later times, at a given point, IW cancel out due to their mutual destructive interference and are superceded by small-amplitude, non-propagating buoyancy oscillations. (Below we discuss the application of the impulsive point theory in the considered case in more detail.)

Figure 7 shows also that the distribution of the density fluctuation in the vicinity of the jet core has a two-layer structure with local maxima of ' at *z* 0.5 , *y* 2 . Now it is well known that the fluid density variations in a pancake vortex arise as a result of the cyclostrophic balance, where the centrifugal force inside the vortex is balanced by a vertical pressure gradient force that is provided by a perturbation of the local density (cf. Beckers et al. (2001)). This balance inside the individual vortex deflects isopycnals toward the vortex center and thus creates a two-layer distribution of the local density perturbation. Figure 6 shows that the centers of pancake vortices are located in the horizontal central plane at *y* 2 . Thus it is plausible to assume that the local maxima of ' at *z* 0.5 , *y* 2 in the distribution of the density fluctuation (cf. Fig. 7, *Nt* = 60, 90) are due to the cyclostrophic balance inside the pancake eddies.

Fig. 7. Contours of the *x*-averaged density fluctuation (left column), and horizontal (centre column) and vertical (right column) velocity fluctuations, ' *U <sup>h</sup>* and ' *U <sup>z</sup>* , obtained at time moments *Nt* = 15, 30, 60, 90 (from top to bottom). The respective maximum values are: 0.06, 0.12, 0.028 (*Nt* = 15); 0.03, 0.1, 0.0085 (*Nt* = 30); 0.016, 0.09, 0.0028 (*Nt* = 60); 0.0095, 0.07, 0.0012 (*Nt* = 90).

vertical. At later times, at a given point, IW cancel out due to their mutual destructive interference and are superceded by small-amplitude, non-propagating buoyancy oscillations. (Below we discuss the application of the impulsive point theory in the

Figure 7 shows also that the distribution of the density fluctuation in the vicinity of the jet core has a two-layer structure with local maxima of ' at *z* 0.5 , *y* 2 . Now it is well known that the fluid density variations in a pancake vortex arise as a result of the cyclostrophic balance, where the centrifugal force inside the vortex is balanced by a vertical pressure gradient force that is provided by a perturbation of the local density (cf. Beckers et al. (2001)). This balance inside the individual vortex deflects isopycnals toward the vortex center and thus creates a two-layer distribution of the local density perturbation. Figure 6 shows that the centers of pancake vortices are located in the horizontal central plane at *y* 2 . Thus it is plausible to assume that the local maxima of ' at *z* 0.5 , *y* 2 in the distribution of the density fluctuation (cf. Fig. 7, *Nt* = 60, 90) are due to the cyclostrophic

U'

<sup>h</sup> U'z

max

0

Y

Fig. 7. Contours of the *x*-averaged density fluctuation (left column), and horizontal (centre column) and vertical (right column) velocity fluctuations, ' *U <sup>h</sup>* and ' *U <sup>z</sup>* , obtained at time moments *Nt* = 15, 30, 60, 90 (from top to bottom). The respective maximum values are: 0.06, 0.12, 0.028 (*Nt* = 15); 0.03, 0.1, 0.0085 (*Nt* = 30); 0.016, 0.09, 0.0028 (*Nt* = 60); 0.0095, 0.07,


considered case in more detail.)

balance inside the pancake eddies.

Z




0.0012 (*Nt* = 90).



Fig. 8. Instantaneous contours of the vertical velocity *U <sup>z</sup>* in the vertical (*x,z*) plane at *y* = 0 (left column) and *y* = 1.5 (right column) at time moments *Nt* = 15 (top frame) and *Nt* = 90 (bottom frame).

Figure 9 presents the distribution of the vertical velocity, *Uz* , in the horizontal (*x, y*)-plane at the distance *z* = 3 above the jet streamwise axis obtained in DNS at time moments *Nt* = 12, 15, 18. The figure shows that the IW patterns in that plane are similar to the IW patterns of stationary lee-waves and characterized by a well-defined streamwise spatial period of about 5 dimensionless units (or about 4-5 sphere diameters in the dimensional form). A similar isophase pattern was observed by BCH in the laboratory experiment in a horizontal plane at a distance of three sphere radii away from the towing axis. Note also that this spatial period is close to the wavelength of the spiral instability mode of the non-stratified jet flow 2 / 4.3 *s s k* for 1.46 *<sup>s</sup> k* (Batchelor and Gill 1962).

Figure 10 shows the temporal development of the vertical velocity and density, *Uz* and *ρ*, obtained in DNS at the point located at *x* = 18, with the angle with respect to the vertical *θ* =

Internal Waves Radiation by a Turbulent Jet Flow in a Stratified Fluid 111

x

Fig. 9. Instantaneous contours of the vertical velocity *Uz* in the horizontal (*x,y*) plane above the jet streamwise axis at *z* = 3 at time moments *Nt* = 12, 15, 18 (from top to bottom).

In the considered case, an IW radiation source function can be evaluated from the equation for the vertical velocity derived by Phillips (1977). In this equation, the non-linear terms are not discarded but can be regarded as an effective source of IW. This equation is derived in

<sup>2</sup> ( , , ,) *U N U Qxyzt z hz <sup>t</sup>*

<sup>2</sup> ( , , ,) *<sup>z</sup> j hj <sup>j</sup>*

*x zt x xt x*

2

full and horizontal Laplasian operators, respectively; summation over the repeated indexes

*U U <sup>b</sup> Qxyzt U U U*

(4.1)

*j jj*

222

2 22 *xyz* 

,

0 5 10 15 20 25 30 35

3

2 22

2


(4.2)

are the

2 2

*<sup>h</sup>* 2 2 *x y* 

2

Uz

y




the inviscid case in the form:

where the IW source function is

In (4.1) and (4.2) the notations are as follows:

0

5

0

5

0

5

<sup>0</sup> 40 at the distance of 2 2 *y z* 3 from the jet streamwise axis for Richardson numbers Ri = 1 and Ri = 3. The figure shows that in both cases IW packet arrives at the observation point at time moment 10 *Nt* which is independent of the Richardson number. At the given point, the IW paket amplitude grows at times 20 *Nt* whereas the velocity *Uz* and density *ρ* oscillate with a phase shift / 2 typical of a monochromatic, linear internal wave (Phillips 1977). At times *Nt* > 30 the oscillations of *Uz* and *ρ* become incoherent and IW packet amplitude decreases. At sufficiently late times (*Nt* > 60) there remain only small-amplitude oscillations of the velocity and density with the buoyancy frequency *N*. These smallamplitude oscillations are similar to the buoyancy oscillations which supercede internal waves, which are mutually cancelled due to the destructive interference as their wavelengths become of the order of the source diameter (i.e. of order one in the considered case) (Voisin 1991).

Figure 11 presents spatially averaged frequency spectra of the vertical velocity obtained in DNS for Richardson numbers Ri = 1 and 3. Each spectrum was obtained by Fourier transform of time series ( ) *U t <sup>z</sup>* and averaged over 100 points located uniformly at a cylindrical surface with radius 2 2 *y z* 3 in the range 0 0 0 90 . The figure shows that the velocity spectra are characterized by a well pronounced peak at the frequency / 0.7 *N* which corresponds to the IW propagation angle <sup>0</sup> arccos / 45 *N* . This value is close to the prediction of the linear theory stating that the amplitude of the waves radiated impulsively by a small-size source has a maximum at <sup>0</sup> arctan 2 55 . This is also in general agreement with the data in Figs. 7 and 8 discussed above.

Figure 12 presents the dependence of the internal waves amplitude ( *iw* ) vs. the inverse Richardson number ( <sup>1</sup> *Ri* ) obtained in DNS for two different values of the Reynolds number, Re = 400 and 700. The IW amplitude was evaluated in each DNS run as a maximum of the r.m.s. density fluctuation '( , , ) *y z t* (2.16) at the distance 2 2 *y z* 3 from the jet streamwise axis. (Note that since we consider the linear reference density profile, '( , , ) *yzt* is also equivalent to the average fluid particle vertical displacement amplitude i.e. the average IW amplitude.). The DNS results show that, for all considered Ri, the maximum of ' is observed at time moment *Nt* ≈ 15 in the range 0 0 40 60 . Figure 12 shows that IW amplitude decreases with increasing Ri. The figure also shows an asymptotic estimate for the amplitude of the impulsively emitted internal waves (in dotted line).

Below in Section 4 the application of the impulsive source theory to explain the observed IW kinematics and dynamics is discussed.
