September 2014 

NASA/TM–2014-218523 

Direct Numerical Simulation of a Weakly 

Stratified Turbulent Wake

J. A. Redford 

University of Southampton, Southampton, United Kingdom 

T. S. Lund 

NorthWest Research Associates, Inc., Boulder, Colorado 

Gary N. Coleman 

Langley Research Center, Hampton, Virginia 
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page 

at http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443-757-5803 

• Phone the NASA STI Information Desk at 

443-757-5802 

• Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

September 2014 

NASA/TM–2014-218523 

Direct Numerical Simulation of Weakly 

Stratified Turbulent Wake 

J. A. Redford 

University of Southampton, Southampton, United Kingdom 

T. S. Lund 

NorthWest Research Associates, Inc., Boulder, Colorado 

Gary N. Coleman 

Langley Research Center, Hampton, Virginia
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

Acknowledgments 

This project was supported by the Revolutionary Computational Aerosciences project under 

NASA’s Fundamental Aerodynamics program, by the UK Ministry of Defence (Dstl project 

RD027-02854, managed by Dr R.P.Hornby) and by the UK Turbulence Consortium, funded by 

the UK Engineering and Physical Sciences Research Council (EPSRC), under Grants 

EP/D044073/1 and EP/G069581/1. It made use of the facilities of HPCx and HECToR, the UK’s 

national high-performance computing service, which is provided by UoE HPCx Ltd at the 

University of Edinburgh, Cray Inc and NAG Ltd, and funded by the Office of Science and 

Technology through EPSRC’s High End Computing Programme. We are grateful to Dr. S. 

L.Woodruff for reviewing the manuscript. 

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 

official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics 

and Space Administration. 
Abstract

Direct numerical simulation (DNS) is used to investigate a time-dependent turbulent wake

evolving in a stably stratified background. A large initial Froude number is chosen to allow

the wake to become fully turbulent and axisymmetric before stratification affects the spreading rate of the mean defect. The uncertainty introduced by the finite sample size associated

with gathering statistics from a simulation of a time-dependent flow is reduced, compared to

earlier simulations of this flow. The DNS reveals the buoyancy-induced changes to the turbulence structure, as well as to the mean-defect history and the terms in the mean-momentum

and turbulence-kinetic-energy budgets, that characterize the various states of this flow –

namely the three-dimensional (essentially unstratified), non-equilibrium (or ‘wake-collapse’)

and quasi-two-dimensional (or ‘two-component’) regimes observed elsewhere for wakes embedded in both weakly and strongly stratified backgrounds. The wake-collapse regime is

not accompanied by transfer (or ‘reconversion’) of the potential energy of the turbulence

to the kinetic energy of the turbulence, implying that this is not an essential feature of

stratified-wake dynamics. The dependence upon Reynolds number of the duration of the

wake-collapse period is demonstrated, and the effect of the details of the initial/near-field

conditions of the wake on its subsequent development is examined.

1 Background and objectives

This study is concerned with the effects of stable buoyancy upon the evolution of a turbulent

wake. These effects can be parameterized by the Froude number F = U/NL, where U and

L are relevant velocity and length scales and N = [−g(dρ/dz)/ρ∞]

1/2 is the Brunt-V ̈ais ̈al ̈a

frequency (where g is the gravitational acceleration, z is the vertical coordinate, ρ is the

mean background density and ρ∞ is a characteristic/reference density). The Froude number

is an inverse measure, being small when buoyancy is important and large when it is not.

For the flow over a wake-generating body, an appropriate Froude number can be defined as

F∞ = 2U∞/ND, where U∞ and D are the speed and effective cross-sectional diameter of the

wake-generating body. A local, time-dependent, alternative is the deficit Froude number,

Fd = 2Ud/N δ, where Ud and δ are the maximum mean wake-velocity deficit and mean wake

width. Since the velocity deficit decreases with distance from the body while the wake width

increases, the wake Froude number strictly decreases with distance. This implies all wakes

in stratified fluids will be affected by buoyancy at some point in their development.

Here we consider the weakly stratified case with initially high Froude number (Fd 1),

such that the wake is largely unaffected by buoyancy for some distance downstream of

the body. Weakly stratified wakes evolve through three distinct regimes, which Spedding

[1] referred to as the three-dimensional (3D), non-equilibrium and quasi-two-dimensional

states. A review of the towed-wake experiments from which this picture derives is given

by Diamessis, Spedding & Domaradzki [2], while a model that recovers the three regimes,

along with two others, is presented by Meunier, Diamessis & Spedding [3].

During the 3D regime the high-F∞ wake is initially unaffected by buoyancy and thus

produces fully three-dimensional turbulent motions for a period of time, while the wakedefect Froude number remains greater than order 1. The mean and turbulence statistics

are consistent with unstratified flows during this time and exhibit self-similar behavior

indistinguishable from that observed in unstratified wakes (Tennekes & Lumley [4], Bevilaqu

& Lykoudis [5]). This situation is in contrast to flows with initial Fd of order one or smaller,

which are immediately altered by stratification, so that the canonical self-similar state is

never realized. In the extreme initial-Fd → 0 case, buoyancy suppresses vertical motion

to the extent that three-dimensional turbulence never appears (Chomaz et al. [6], Riley &

Lelong [7]). For towed-sphere experiments, the critical value below which the 3D regime

1
does not exist is F∞ ≈ 4; cf. Dommermuth et al., [8].

When the initially three-dimensional turbulence in the high-F∞ wake begins to ‘feel’ the

effects of buoyancy, the non-equilibrium (NEQ) (also known as ‘vertical collapse’) regime

begins. (The term ‘non-equilibrium’ can be viewed as a synonym for ‘intermediate’ or

‘transitional’, and should not be interpreted literally: it does not imply, for example, that

in the pre-NEQ state the rates of production and dissipation of turbulence kinetic energy

balance, since this type of equilibrium does not happen even in the purely unstratified selfsimilar benchmark.) When the wake enters the NEQ phase, it ceases to spread vertically;

in fact, the height has been observed in some cases to slowly decrease with time during

this period (Brucker & Sarkar [9]; Diamessis et al. [2]; Lin & Pao [10], the latter for the

self-propelled case), although this may be at least in part a consequence of how the height

is defined. As we will see below, a distinction should be made between the heights based on

the mean defect (typically measured using a Gaussian curve fit) and on the extent of the

wake turbulence (revealed for example by shadowgraph or dye visualizations).

As the wake stops spreading vertically, its lateral growth continues at roughly the same

rate observed during the 3D phase (Diamessis et al. [2]). (We reserve the term ‘lateral’ to

indicate the direction aligned with the axis orthogonal to the streamwise-vertical plane.)

Whether the rate of lateral spreading remains exactly or only approximately constant as

the flow passes from the 3D to NEQ states is an open question. Brucker & Sarkar [9] cite

different power-law exponents for the mean-defect width during the different regimes, while

Diamessis et al. [2] find t

1/3 behavior for the entire wake history.

The NEQ period is also characterized by a large reduction in the rate of decay |dUd/dt|

of the maximum mean defect velocity Ud, and also the rate of decay of the turbulence

kinetic energy. Spedding [1] and others (e.g. Dommermuth et al. [8]; Diamessis et al. [2])

have described these as ‘re-stratification’ effects, attributed to the transfer of potential

to kinetic energy, as the wake turbulence collapses, and heavy/cold fluid parcels fall, and

light/warm ones rise, towards their original positions. Although this explanation is plausible

and consistent with the behavior of the mean defect, it is not the only possibility, since the

mean and turbulence kinetic energy can also be enhanced by alterations to the production

and dissipation processes (cf. Brucker & Sarkar [9]). We shall present evidence for a picture

of NEQ dynamics that is not based on the net transfer of potential to kinetic energy, which

for the turbulence is consistently in the opposite direction and for the mean does not act

directly upon the streamwise component.

The NEQ phase is observed to begin at roughly two buoyancy periods (N t ≈ 2) for a

wide range of F∞ > 4 (Spedding [1]). Its duration, on the other hand, has been predicted by

Meunier et al. [3], and demonstrated by the DNS of Brucker & Sarkar [9] and the large-eddy

simulations (LES) of Diamessis et al. [2], to increase with Reynolds number R∞ = U∞D/ν.

While the evidence for this Reynolds-number dependence is strong, it remains at present

rather qualitative, since it is influenced by the statistical uncertainty associated with finitedomain numerical simulations (see below).

After the NEQ phase, the wake enters the quasi-two-dimensional (Q2D) state, characterized by the emergence of large-scale coherent structures. (Godeford & Cambon [11]

point out that ‘quasi-two-component’ is a more accurate description of flows like this, whose

vertical motions have been strongly attenuated by stable buoyancy, since significant vertical

variations of the other two components are not necessarily negligible. Nevertheless, to be

consistent with previous studies, we shall continue to use ‘Q2D’ to refer to this aspect of

the wake evolution.) The towed-sphere measurements of Chomaz et al. [6] revealed that

flattened ‘pancake-like’ eddies eventually appear for both weak and strong stratification (cf.

Spedding [12]). The maximum mean velocity begins to decrease more quickly than it did

during the NEQ regime, while the mean defect spreads laterally at a rate comparable to

the 3D value. The vertical spreading also begins again, growing at a rate consistent with

viscous diffusion (Meunier et al. [3], Brucker & Sarkar [9], Diamessis et al. [2]).

2
While the effects of buoyancy that lead from the 3D to the NEQ and then the Q2D

regimes are well-understood on a qualitative level, there is lesser understanding of the detailed turbulent processes that drive the regime changes. Our aim in the present study is

to use numerical simulation to quantify the mechanisms involved, and thereby provide answers to questions such as the one mentioned above regarding the hypothesized transfer of

potential to kinetic energy during the NEQ regime. Nearly all wake simulations performed

to date have made use of the temporal idealization where a fixed streamwise segment of the

wake is evolved in a reference frame attached to the undisturbed fluid. The weak streamwise

variation of statistics is ignored so that periodic boundary conditions can be used in the

streamwise direction. This configuration allows for very efficient numerical procedures and

has proven to yield results that are in reasonable agreement with laboratory measurements

(cf. Gourlay et al. [13]). Although we employ the same temporal idealization as those before

us, we have chosen a longer streamwise domain than has been used in prior studies. Since

the streamwise direction is assumed to be homogeneous, it forms the basis for computing

the (time-dependent) turbulence statistics. While relatively short domains require fewer

mesh points and proportionally less computer time, they provide a smaller averaging base

at each time; a domain of minimal size may contain only one or two of the large-scale turbulent structures present at that instant in time. While such a simulation progresses without

difficulty, there can be large uncertainty, which increased in time, in the computed statistics

due to inadequate sample. Here we use a streamwise domain size that is between 1.2 and

5.6 times longer than those employed in previous wake simulations (cf. Gourlay et al. [13],

Brucker & Sarkar [9], and Diamessis et al. [2]). This extended streamwise domain allows

for a greater population of ‘large-eddy samples’ and provides less uncertain statistics. This

strategy leads to resolution of earlier discrepancies, such as that in the lateral spreading rate

during the NEQ regime reported by [9] and [2].

As with the discrepancy in the lateral wake spreading rate, there is no consensus on

the Reynolds-number dependence of the duration of the NEQ regime. This feature is also

influenced by statistical uncertainty, and prior results reported by Brucker & Sarkar [9] and

Diamessis et al. [2] have relatively high uncertainty due to limited large-eddy samples. The

increased domain size and a separate numerical experiment (discussed in detail below) will

be used to demonstrate the factors affecting the time at which the NEQ state begins to

transition to the Q2D regime.

Finally, the initial-condition dependence of the wake development will be examined, in

light of the non-universality of the self-similar states observed for late axisymmetric wakes

in unstratified media. The non-universality of the unstratified case is exemplified by mean

spreading rates and turbulence-to-mean ratios that vary with the shape of the wake-creating

object (Bevilaqua & Lykoudis [5]) (although they eventually, at distances corresponding to

1000s of body sizes downstream, relax to a universal value; see Redford et al. [14], referred

to hereinafter as RCC). This non-universality has been traced to differences in the near-field

turbulence structures shed from bodies of different shape. For the stratified counterpart,

the available data suggest that the structures of the early and late wake are unrelated, with

pancake eddies ultimately resulting from a wide variety of initial states, some of which are

non-turbulent or even non-wake-like [2,13,15,16]. Furthermore, Spedding [1] showed that for

sphere wakes, nominal universality is obtained for the mean-defect history Ud(t)/V∞ when

weighted by F2/3 ∞ . This scaling was generalized to bodies of various shapes by Meunier &

Spedding [17] (see also [3]) by accounting, via an effective diameter, for the net momentum

extracted by the body. The founding assumption of their model will be assessed below.

The issues raised above will be addressed using DNS of a wake at high F∞ and relatively high R∞. To put this simulation in context, we mention the earlier numerical

studies of stably stratified wakes generated by constant-velocity axisymmetric towed bodies

performed, for high-F∞, by Gourlay et al. [13] (via DNS), Diamessis et al. [2] (LES) and

3
Brucker & Sarkar [9] (DNS), the latter as part of a study of differences between wakes downstream of towed- and self-propelled bodies in high- and low-F∞ backgrounds (de Stadler

& Sarkar’s [18] DNS includes the effect of body acceleration for moderate stratification).

Lower F∞ cases have been simulated by Dommermuth et al. [8], by LES, and (for a range

of Prandtl numbers) de Stadler et al. [19], by DNS. Provided it can reach relevant Reynolds

numbers (and produce accurate statistics), a numerical approach has many advantages, supplying detail not available from experiment. For example, Gourlay et al.’s , [13] DNS revealed

the vortex-ring structure underlying the Q2D pancake eddies, while the LES of Diamessis

et al. [2] documented the presence of small-scale Kelvin-Helmholtz instabilities during the

NEQ regime and its Reynolds-number dependence. Dommermuth et al.’s [8] low-F∞ LES

quantified the relative importance of the terms in the mean- and turbulence-kinetic-energy

budgets, during the NEQ and Q2D regimes, that ensue when the turbulence is immediately affected by buoyancy. Similar exercises were performed by Brucker & Sarkar [9] and

deStadler & Sarkar [18] for their towed- and self-propelled-body flows.

To summarize, we will use DNS to:

1. Obtain results with improved ‘eddy sample size’ (via increased streamwise domain

length), to reduce the uncertainty of the statistics, especially the mean wake spreading

rate.

2. Assess more completely how the effects of buoyancy alter the wake’s mean momentum

and turbulence energy budgets.

3. Determine whether or not there is a net transfer of potential to kinetic energy during

the wake collapse in the NEQ regime.

4. Quantify the Reynolds-number dependence of the duration of the NEQ regime.

5. Demonstrate how the wake development is affected by the details of the initial condition used in the numerical simulation, comparing with earlier simulation results to

explore the extent to which mean-defect histories exhibit universality.

In § 2, the numerical approach, initialization strategy, and DNS parameters are described.

Section 3 documents the wake development in terms of flow visualisations, mean-flow statistics, and the budgets of mean momentum, and mean and turbulence kinetic energy. Closing

comments and conclusions are summarized in § 4.

2 Problem formulation, numerical strategy and parameters

Following Gourlay et al. [13] (cf. [2, 8, 9, 18, 19]), we (as in RCC) invoke a Galilean transformation, and emulate the actual, spatially developing wake with a parallel-flow, temporally

developing idealization. The resulting streamwise homogeneity allows efficient DNS with a

standard triply periodic FFT-based pseudo-spectral method [20].

For the high-Froude-number flow studied here, vertical density variations are negligible compared to the background reference density ρ∞, allowing stable stratification to be

accounted for with the Boussinesq (divergence-free velocity) approximation. Molecular diffusion of the (temperature-dependent) density field is matched to that for the velocity field

by assuming the Prandtl number ν/α = 1 (cf. [19]), where ν is the kinematic viscosity and

α is the thermal diffusivity, both of which are assumed constant in space and time (and

therefore independent of temperature). The pressure variable in the incompressible NavierStokes equations is eliminated using a velocity-vorticity formulation [21]. A low-storage

4
third-order temporal discretization [22] is applied to the nonlinear terms, while the viscous

terms are treated analytically [23].

The coordinate x = x1 denotes the streamwise/longitudinal direction, along the wake

axis, with u = u1 the corresponding streamwise velocity, defined such that positive values

point ‘upstream’ and the maximum occurs at the centerline (figure 9). Because of the

parallel-flow assumption, the uniform background/freestream velocity V∞ amounts to an

arbitrary Galilean shift; we chose V∞ = 0, which implies u → 0 far from the wake centerline.

The lateral (‘spanwise’) and vertical directions are respectively y = x2 and z = x3 (positive

upwards), with velocities v = u2 and w = u3, such that x-y planes are normal to the gravity

vector gi = −gδi3.

Before discussing the numerical parameters, we describe the three types of averaging used

in this paper. The Reynolds decomposition is defined with respect to the mean (denoted by

an overbar) formed by averaging over the homogeneous streamwise direction x; this mean

thus varies with lateral and vertical position, and time: e.g. u = u(y, z, t). (Note that,

unlike for the unstratified wake simulated in RCC, averaging over the azimuthal direction is

not appropriate.) A prime (

) indicates a deviation from this mean, such that the Reynolds

stresses, density-fluctuation variance and turbulent buoyancy flux are written respectively

as −u

iuj, ρ

ρ and −ρw. As in Spedding [1], we also examine histories of bulk averages

of Reynolds-averaged quantities, taken across the inhomogeneous directions (y and z) of

the mean defect, over the region where u > 0.05 Ud (Spedding’s mean was over horizontal

and vertical planes where u > 0.2 Ud). These averages are indicated with angle brackets:

k = (1/2)u

iui = (1/2) Λz

0

Λy

0 ψ u

iui dy dz/Ψ, where ψ = ψ(y, z)|

u/Ud>0.05 is unity if

u(y, z) > 0.05 Ud and zero otherwise, and Ψ = Λz

0

Λy

0 ψ dy dz. A third type of statistic

(indicated by braces { }) is formed from cross-wake integration (not averaging) over the

entire y-z domain, of Reynolds-averaged quantities; examples include the total/integrated

turbulence kinetic energy (TKE) {k} = (1/2) Λz

0

Λz

0 u

iui dy dz and its rate of production

{Ptke} = − Λz

0

Λz

0 (u

iuj∂ ui/∂xj ) dy dz.

Based on our earlier experiences with DNS of unstratified turbulent wakes (RCC), the

wake is initialised as a sequence of vortex rings, identical to Case VR of RCC apart from

the ring spacing (see their figure 1a), embedded in a stably stratified background with

imposed constant/uniform mean-density gradient Γ = (dρ/dz)∞ < 0. This strategy was

found to be an efficient way to produce fully developed wake turbulence containing mature

turbulent structures, which – while not necessarily corresponding to those generated by a

realistic wake-creating body in a real-world application – allows us to perform a numerical

experiment that can be used to address the issues raised above (such as the effect of the

details of the initial conditions on late-wake development).

In order to simulate the wake using periodic boundary conditions, the transported density

variable (= ρ(z, t)+ρ

(x, y, z, t)) is defined as the perturbation from the background density

field ρ∞ + zΓ (where the origin of the z-axis is, say, in the center of the computational

domain).

As in Shariff [24], each ring in the sequence is represented as a circular loop, of radius r0

in y-z planes, of Gaussian azimuthal core vorticity ωθ, with ωθ = (γ/πδ2

ω) exp 

−(s/δω)2



,

where γ is the ring’s circulation, s2 = (x − x0)2 + (r − r0)2, r2 = (y − y0)2 + (z − z0)2 and

(x0, y0, z0) is the location of the centre of the ring in question. The length-scale δω defines

the thickness of the vortex core (chosen as δω = 0.4r0). The rings are spaced uniformly every

2πr0/3 along the wake axis (x, y0, z0) (RCC used 16πr0/27 ≈ 0.59πr0); the rings combine

to induce an initial mean velocity defect of maximum velocity Udo and width ho, the latter

defined as the radial distance from the wake centreline to the location at which the mean

streamwise velocity u is half of Ud.

The wake turbulence is the result of a nonlinear/bypass transition, which is triggered

by a Δr = ±0.004r0 maximum radial displacement of the ring circumference, in the y-z

5
plane, for each of the azimuthal modes 5 − 32, with random azimuthal phases, using different random-number seeds for each ring. As the individual vortex rings break down and

disturb their neighbours, a vericose-type instability ensues (see figure 5a of RCC), leading

to fully developed turbulence relatively quickly – much sooner than, for example, when

very-small-amplitude disturbances are added to a mean defect profile, such that transition

is via growth of the helical linear-instability modes (as was done for Case SD of RCC). The

approach used by Gourlay et al. [13] also avoids the need to simulate the natural linearinstability transition process, by prescribing a Gaussian mean velocity and superimposing

random finite-amplitude perturbations whose magnitudes were set by matching a desired energy spectrum, and then clipped/filtered in physical space with an axisymmetric/cylindrical

window. Similar procedures were followed by Brucker & Sarkar [9] and Diamessis et al. [2],

although the former’s DNS involved an initial adjustment period during which the mean

flow was held fixed as the fluctuations developed, and the latter’s LES required a precursor

‘relaxation’ simulation, during which the effective gravity was slowly increased. Despite

being necessarily somewhat ‘synthetic’ in nature, those initialization schemes all produced

statistics that very quickly exhibit expected behavior (see figure 21). The primary difference

between the present and earlier strategies is that here coherent structures are ‘hard wired’

into the initial-condition turbulence.

The initial Froude number specified here is quite large, to ensure the wake turbulence

is fully developed and physical when it is first affected by the stratification. The initial

mean-defect Froude number Fdo = Udo/Nho = 23.7; this corresponds to F∞ = 2V∞/ND of

approximately 130, assuming the uniform flow speed V∞ ≈ 5Udo and virtual-body diameter

D ≈ 2ho (see RCC). The initial mean-defect Reynolds number Rdo = Udoho/ν = 2625,

such that the equivalent R∞ = V∞D/ν ≈ 2.6 × 104. This is well above the critical value

of R∞ ≈ 5 × 103, found by Spedding, Browand & Fincham [25] to be required for the

turbulent wake structures and decay rates to exhibit self similarity – although it is not

above the threshold (if one exists) beyond which all features (such as the duration of the

various self-similar regimes) can be expected to be independent of Reynolds number (cf.

Diamessis et al. [2]).

A modified version of this case will also be considered, which is equivalent except that

beginning from N t = 1.9 the viscosity varies at each time step, such that the mean-defect

Reynolds number, based on the integral scale r∗ = ((ln 2/π)(Id/Ud))1/2, where Id = {u}

is the volume-flux deficit, remains fixed thereafter at Ud r∗/ν = 1500; see + symbol in

figure 1c, d. (The regridding/domain-expansion strategy described below was also employed,

although, towards the end of the run, not at identical times.) Results from this ‘artificial’

numerical experiment (cf. [26]), which we shall refer to as Case R1500, will be compared

to those from the primary (constant-viscosity/variable-Renolds-number) DNS, to determine

the Reynolds-number dependence of the mean-defect histories; unless otherwise stated (see

figure 7) all the new simulation results shown below will be from the primary DNS. Two

additional test runs were also performed, to explore numerical-validation issues, as explained

below.

For Gourlay et al. [13], Fdo = 34.0 and Rdo = 1472. The more recent simulations tended

to use lower F and higher R: Brucker & Sarkar’s [9] Case TR50F04 (i.e. F∞/2 = V∞/ND =

4, R∞ = 50 000) corresponds to Fdo ≈ 0.75 and Rdo ≈ 3250 (they assume Udo/V∞ = 0.11

and D ≈ 1.7ho). They also simulated a lower-R, higher-F flow with which we will compare

(their ‘TR10F20’), with (Fdo, Rdo) ≈ (3.75, 650). Comparison will also be made with the

recent F∞ = 64 (Fdo ≈ 14.5), R∞ = 105 (Rdo ≈ 1.1 × 104) LES of Diamessis et al. [2], for

which the eddy-resolving capability was comparable to the present DNS. In fact, although

the smallest resolved scales are not those that actually appear at R∞ = 105, the LES has

the advantage of taking full benefit of the available spatial resolution – since the magnitude

of the subgrid-scale viscosity of the LES can be smaller than the molecular viscosity of the

DNS – and thus has the capacity to capture scales smaller than those of an over-resolved

6
    











Ud/Nhz

(a)

    











Ud/Nhy

(b)

    







N t

Ud hz/ν

(c)

    







N t

Ud hy/ν

(d)

Figure 1. Mean (a, b) Froude number Fd = Ud/Nh and (c, d) Reynolds number Rd =

Udh/ν histories: , present DNS [(Fdo, Rdo) = (23.7, 2223)]; +, time at which the

variable-viscosity run (Case R1500) first reaches its minimum Reynolds number, Ud r∗/ν =

1500 (see discussion regarding figure 7); ◦, Gourlay et al. [13] [(Fdo, Rdo) ≈ (34, 1472)]; ,

Case TR50F04 of Brucker & Sarkar [9] [(Fdo, Rdo) ≈ (0.75, 3250)]; •, Case TR10F20 of

Brucker & Sarkar [9] [(Fdo, Rdo) ≈ (3.75, 650)]; , Case R100F64 LES of Diamessis et al. [2]

[(Fdo, Rdo) ≈ (14.5, 11 000)].

7
DNS at smaller Reynolds number, such as is used here.

The variations of the mean-defect-based Froude and Reynolds numbers for the present

and earlier simulations are shown in figure 1, defined in terms of both the wake’s mean height

hz and width hy. Evidence of each of the 3D, NEQ and Q2D regimes can be seen. The

tendency for the Froude numbers to all fall monotonically indicates that all the flows, even

those with very large Fdo, eventually are dominated by buoyancy effects. The anisotropy

introduced by the stable stratification can be inferred in the differences between the Froude

and Reynolds numbers based on hz and hy. The wake evolution is such that the power-law

decay of the height-based Reynolds number Ud hz/ν varies little over the life of the flow,

while Ud hy/ν remains relatively constant, especially during the NEQ phase (cf. Diamessis et

al. [2]). The time required (t ≈ 1/N) for the vortex-ring initialization (solid line) to produce

a fully turbulent state is apparent, which can be contrasted with the approaches used to

define the N t = 0 state for the other simulations (recall that all involved specifying finite

disturbances typical of the fully developed wake downstream of a sphere). The question of

appropriate virtual origins, needed when comparing the present and earlier cases, will be

taken up below.

To prepare for the possibility of the wake generating significant internal/buoyancy waves,

a ‘fringe’ region is adopted at the lateral and vertical edges of the periodic domain, to prevent

the energy radiated in the waves from re-entering the domain, and spuriously affecting the

wake dynamics. A forcing term F (similar to that used by de Stadler et al. [19] and

de Stadler & Sarkar [18]) is added to the right-hand-side R of the momentum and scalar

equations

∂q

∂t = R + F, (1)

which, in the fringe zone, drives q = (u, v, w, ρ

)T to a mean reference qˆ, defined as its

average along the edge of the domain. The fringe forcing is written as F = −σ(y, z)(q − qˆ)

where the fringe function

σ(y, z) = Cσ [ σy(y) + σz(z) − σy(y)σz(z)] , (2)

for 0 ≤ y ≤ Λy and 0 ≤ z ≤ Λz, with Cσ = 0.04 and

σy(y)=1 − 1

2



tanh af

2y

λy

− 1

	 + tanh 

af

2(Λy − y)

λy

− 1

	
 , (3)

σz(z)=1 − 1

2



tanh af

2z

λz

− 1

	 + tanh 

af

2(Λz − z)

λz

− 1

	
 , (4)

where af = 2.5, λy = 30Λy/Ny, λz = 30Λz/Nz and Λy and Λz are respectively the lateral

and vertical domain sizes.

The effectiveness of this scheme was tested by applying it to the internal wave system

created by a point vortex (a multiple-frequency analogue of the vibrating-cylinder, ‘Saint

Andrew’s cross’, experiments of Mowbray & Rarity [27] (also in [28] and [29]), to verify that

the internal waves neither re-entered the domain or reflected off the fringe boundaries. The

energy absorbed by the fringe was calculated during the turbulent wake simulation. The

strongest activity was detected during the NEQ phase of the simulation, when animations of

quantities such as ∂w/∂z (cf. figure 12b) suggest that the waves are generated by individual

turbulent structures (rather than, say, mean collapse of the wake defect associated with entry

into the Q2D regime). For the present high-F flow, the internal-wave energy was found to

be only a small percentage of the overall dissipation (see discussion regarding figure 3b).

Another computational challenge associated with this flow is defining an appropriate

domain. The streamwise two-point correlations in figure 2a indicate that the streamwise

domain size Λx is adequate throughout the simulation, and individual structures fit comfortably within the domain (since their streamwise integral scale is much smaller than Λx),

8
     	

	



	



rx/Λx



R11

(a)

    	










	



κx



E11

(b)

Figure 2. Evolution of (a) streamwise two-point autocorrelations and (b) streamwise spectra

of streamwise velocity: /◦, N t = 0.2 (in b only); /, N t = 3.2; /,

N t = 21; /•, N t = 136; /, N t = 587. Straight solid line in (b) indicates

κ−5/3 x behavior. Spectra E11 and correlations R11 = Q11/Q11(0) are from averages of

E11 and Q11 over lateral and vertical planes though the wake centerline in regions for

which u ≥ 0.5 Ud, where u

u = 

κx E11(κx) and Q11(rx) = u

(x)u(x + rx). Spectra and

wavenumbers in (b) normalized by Ur = 2.09Udo and Lr = ho.

even at very late times. The large-scale waviness during the late stages is due to the relatively low Reynolds number (i.e. narrow range of streamwise scales), and especially the

periodicity of the large pancake-eddy structures present in the Q2D regime (figure 6). However, the streamwise domain size also directly affects the quality of the mean statistics (and

possibly the dynamics of the merging of the large-scale structures), since it defines the size

of the ‘large-eddy sample’ (which will decrease as time passes) that enters the streamwise

mean. Therefore, the choice was made to have a relatively large streamwise extent, with

Λx/ho = 96π ≈ 302 (cf. Gourlay et al.’s [13] Λx/ho = 245, Brucker & Searcher’s [9]

Λx/ho ≈ 105, and Diamessis et al.’s [2] Λx/ho ≈ 54).

The possibility that the height and especially the width of the wake will outgrow the

periodic domain is of even greater concern. To avoid this problem, without unduly wasting

computational resources, the lateral and vertical domain sizes Λy and Λz are dynamically

increased at appropriate times during the run, by projecting the periodic solutions onto a

domain of greater extent in y, and lesser in z, with the same grid spacing. Details of this

re-projection procedure are given in RCC. The spatial resolution is also reduced at certain

times during the simulations, as allowed by the monotonically decreasing Reynolds number

(see table 1).

The choices for the lateral and vertical domain are tabulated in table 1, with respect to

both the initial ho and local mean lateral hy and vertical hz half-widths. The criterion found

by RCC for the unstratified case, that the lateral and vertical width h must remain less than

0.17 of the corresponding domain size, is satisfied throughout the simulation. Confirmation

that stratification does not significantly alter this guideline is provided in figure 3a, where

the DNS histories are given along with those from a test run identical to the case defined

in table 1, except that the lateral and vertical domain sizes remain fixed at their values at

N t = 28 (corresponding to the end of the NEQ regime). There is good agreement between

the two runs (including for the maximum Ud, width hy and height hz of the mean defect,

as well as the streamwise u

, lateral v and vertical w velocity fluctuations) until N t ≈ 340,

at which point hy/Λy ≈ 0.15 and hz/Λz ≈ 0.12 for the smaller/fixed-domain results. These

9
Nt Mx Δx/ηmin My Δy/ηmin Λy/ho hy/Λy Δy/hy Mz Δz/ηmin Λz/ho hz/Λz Δz/hz

0 3072 3.8 256 1.0 2π 0.159 0.02 256 1.0 2π 0.159 0.02

0.87 7.0 1.7 0.155 0.02 1.7 0.157 0.020.87 384 3π 0.103 0.02 384 3π 0.105 0.02

1.12 6.7 256 2.5 0.112 0.02 256 2.5 0.110 0.021.12 384 4.5π 0.075 0.02 384 4.5π 0.073 0.02

2.86 2.9 256 1.6 0.131 0.02 256 1.6 0.120 0.022.86 384 6.75π 0.087 0.02 384 6.75π 0.080 0.02

5.38 1.8 256 1.5 0.132 0.02 256 1.5 0.087 0.035.38 384 10.13π 0.088 0.02 0.087 0.03

9.72 2560 1.4 256 1.5 0.108 0.02 192 1.3 0.085 0.049.72 384 15.19π 0.072 0.02 0.085 0.04

18.6 2304 1.0 256 1.4 0.078 0.03 0.8 0.078 0.0518.6 384 22.78π 0.052 0.03 0.078 0.05

28.1 1536 1.2 320 1.4 0.057 0.04 128 1.0 0.077 0.0728.1 480 34.17π 0.038 0.04 192 10.13π 0.052 0.07

35.4 1024 1.7 384 1.6 0.040 0.04 1.0 0.051 0.0746.1 512 3.3 256 2.4 0.045 0.06 128 1.4 0.051 0.1066.0 384 4.3 192 3.0 0.051 0.07 96 1.8 0.052 0.1366.0 256 45.56π 0.038 0.07 128 13.5π 0.039 0.13

109 3.5 192 3.4 0.066 0.05 96 2.0 0.049 0.14109 256 60.75π 0.050 0.05 128 18.0π 0.037 0.14

215 256 3.1 192 2.6 0.072 0.05 96 1.5 0.056 0.12215 256 81π 0.054 0.05 128 24π 0.042 0.12

587 1.5 1.3 0.081 0.03 0.8 0.067 0.08841 1.3 1.1 0.094 0.03 0.6 0.079 0.07

Table 1. Numerical parameters. The numbers of streamwise Mx, lateral My and vertical Mz Fourier-expansion coefficients are two-thirds of the

corresponding collocation/quadrature points: Nx = 3Mx/2, Ny = 3My/2 and Nz = 3Mz/2. Streamwise Δx, lateral Δy and vertical Δz grid

spacings given by Λx/Nx, Λy/Ny and Λz/Nz, respectively. Minimum Kolmogorov length scale ηmin = (ν3/max)1/4, where max is the maximum

of the streamwise-mean rate of dissipation of TKE at each time. The length Λx of the streamwise domain is fixed at 96πho, where ho is the

mean wake half-width/height at N t = 0, while Λy and Λz increase at times shown (indicated by vertical gaps between rows). Values remain

unchanged until a new entry is given.

10
	



	



    



	











N t

h/ho u/Udo

(a)    







  

	

	



  



	



N t

d{k}/dt (×10−2) (×10−4) (×10−6)

(b)

Figure 3. (a) Comparison of results from the present DNS and a fixed-domain test run

beginning from N t = 28.1, with Λy and Λz maintained at 34.17πho and 10.13πho, respectively, such that the domain is not expanded at N t = 66, 109 and 215; see table 1): lines,

present DNS; symbols, fixed-domain test run (◦, hy; •, hz; , Ud; , u

u1/2; , vv1/2;

, (ww1/2, where streamwise u

u, lateral vv and vertical ww Reynolds stresses

are cross-wake averages over the region in the y-z plane where u > 0.05Ud). (b) Histories

from the present DNS of the rate of change of integrated TKE {k} and net production

minus dissipation plus buoyancy {Ptke}−{tke} + {Btke} (see § 3.5 for definitions): ,

d{k}/dt; , {Ptke}−{tke}; ◦, {Ptke}−{tke} + {Btke}. Budget terms in (b) given in units of

U3

r Lr, where reference velocity Ur = 2.09Udo and length Lr = ho.

ratios are larger than those produced by the expansions used for the NEQ and Q2D phases

of the present DNS, which implies the results have not been adversely affected by a too-small

domain.

Table 1 also shows the changes in resolution over the course of the simulation. The

streamwise Δx, lateral Δy and vertical Δz grid spacings are small enough multiples of the

minimum Kolmogrov scale ηmin throughout the simulation (between 1 and 4, except for a

brief period at very early times when Δx reaches 7ηmin) that we can safely assume the results

faithfully capture all the significant spatial scales of the turbulence. The temporal scales

are resolved by ensuring the Courant-Friedrichs-Lewy (CFL) number does not exceed √2.8.

These conclusions are supported by the streamwise spectra in figure 2b. The signature of

the vortex-ring characterize is apparent at N t = 0.2. The spectra also indicate the moderate

Reynolds number range considered here, in that κ−5/3 behavior is apparent, but only at

early times, bracketing N t = 3. A more empirical check of the resolution is provided by

the histories in figure 3b, comparing the sum of the production, dissipation and buoyancy

terms (the fringe contribution is not included) in the integrated TKE equation (the open

circles) with its integrated rate of change (solid line). The good agreement between the two

points to the accuracy of the present DNS. It also implies that the energy of the internal

waves generated by the turbulence during the NEQ phase is negligible, since the integral of

the flux-divergence terms is essentially zero compared to that of the production, dissipation

and buoyancy. This is in contrast to the low-F case, where departure from the 3D regime

is accompanied by significant TKE transfer away from the wake (Brucker & Sarkar [9]).

We note that figure 3b also reveals the integrated buoyancy term (the difference between

the solid diamonds and open circles) is negligible except between N t ≈ 1.5 and 15, and that

during this time it acts solely as a sink of TKE. This provides an answer to the question

raised in the introduction, regarding the defining NEQ mechanism, since it indicates that

the transition from the 3D to the NEQ regime is not accompanied by a net transfer of

11
potential to kinetic energy.

A final comment regarding the present simulation strategy is warranted: the lateral- and

vertical-domain expansion is accompanied by a coarsening of the resolution across the mean

wake profiles u(y, z0) and u(y0, z), quantified in table 1 by a tendency for Δy/hy and Δz/hz

to increase with time. A check was made that the loss of resolution does not compromise the

shape of the u profiles, by repeating the portion corresponding to the largest Δz (> 0.1hz,

approximately) with double Mz, and confirming that the u results have negligible differences

compared to those shown below.

3 Results

3.1 Flow visualisation

The turbulence structures associated with the 3D, NEQ and Q2D regimes are shown, respectively, in figures 4, 5 and 6. The weak stratification used here allows the wake turbulence to

develop a mature three-dimensional state before buoyancy begins to alter the flow. At the

time shown in figure 4 (N t ≈ 3) the stratification causes the height to be slightly smaller

than the width (see also upper right-hand-side inset plot in figure 21b). However, it is evident that there is a substantial range of scales (cf. figure 2b), and that the wake is to a good

approximation still axisymmetric.

The large-scale lateral oscillation, as well as the flattening and vertical striations, observed in the NEQ structure (figure 5, for N t ≈ 21) are typical of those found in other

experiments and simulations [2, 12, 30], some of which correspond to lower F∞ and higher

R∞. In contrast to the more isotropic eddies seen at N t = 3, here the turbulence in figure 5

contains ‘arrowhead’ structures in the vertical plane through the wake centerline that point

in the direction of motion of the virtual body associated with the wake (cf. figures 10–12 of

Diamessis et al. [2]). A visual estimate of the height of the turbulent region (as opposed to

that of the mean profile) indicates it is similar to, perhaps slightly less than, that for the 3D

regime (compare figure 4b). The change in structure corresponds to enhanced production

of TKE and profound anisotropy in its rate of dissipation (figure 16b). The lateral vorticity

contours in figure 5b suggest the presence of Kelvin-Helmholtz (K-H) billows at the top

and bottom of the wake, in the inclined shear layers of the arrowhead structures (e.g. near

x = 90ho). Diamessis et al. [2] discuss the manner in which the characteristics of these

inclined shear layers depend upon Reynolds number, and the associated question of the how

the layers’ overturning via K-H instabilities can contribute to the generation of turbulence

in stably stratified environments at high Reynolds number. For the present case, we show

below that at the time corresponding to figure 5, viscous effects are not insignificant.

The Q2D regime contains flattened, large-scale, pancake- type motions in horizontal

planes, such as those shown for N t ≈ 587 in figure 6. There are three separate vortical structures that have the two-layered feature that emerged earlier in the simulation (figure 5c),

which bring to mind Godeford & Cambon’s [11] distinction between two-dimensional and

two-component turbulence (see also Spedding [12]). These two-layer eddies are similar

to those observed by Spedding [1] (for N t ≥ 50 and F∞ > 10), by Gourlay et al. [13]

(Nt > 100, F∞ > 150), by Brucker & Sarkar [9] (N t = 125, F∞ = 8), and by Diamessis et

al. [2] (N t = 70, F∞ = 4). Figure 6 is also consistent with the vertical layering of lateral

vorticity fluctuations seen by Spedding [12] (at N t = 64 for F∞ = 4). In light of the vortexring initialization used here, which injects structures very different from the helical modes

associated with the towed-sphere wake (cf. figure 5a of RCC), the present results add to the

growing body of evidence that the Q2D eddies are unrelated to the type of instability by

which the wake became turbulent (see also [2, 13]).

12
(y − y0)/ho

(a)

x/ho

(z − z0)/ho

(b)

Figure 4. Instantaneous flow at N t = 3.2 (Nt

 ̃ = N t − N tvo = 2.5, end of 3D regime;

virtual origin tvo defined in figure 7). (a) Top and (b) side views of isosurfaces of vorticity

magnitude (|ω| = 0.38 Ur/Lr). Mean-flow direction defined as if a virtual wake-creating

body has moved from left to right through the sub-domain shown. Views shown are 8% of

streamwise, and 47% of lateral and vertical, domain.

13
(y − y0)/ho

(a)

(z − z0)/ho

(b)

x/ho

(z − z0)/ho

(c)

Figure 5. Instantaneous flow at N t = 20.9 (Nt

 ̃ = N t − N tvo = 20.2, NEQ regime; virtual

origin tvo defined in figure 7). (a) Top and (c) side views of isosurfaces of vorticity magnitude

(|ω| = 9.7 × 10−4 Ur/Lr). (b) Contours of lateral vorticity ωy in vertical plane through

centerline, in units of 102Ur/Lr. Motion of virtual wake-creating body is from left to right.

Views shown represent subregion of full domain (table 1).

14
(y − y0)/ho

(a)

(y − y0)/ho

(b)

x/ho

(z − z0)/ho

(c)

Figure 6. Instantaneous flow at N t = 586.9 (Nt

 ̃ = N t − N tvo = 586.2, Q2D regime;

virtual origin tvo defined in figure 7). (a) Contours of vertical vorticity ωz in horizontal

plane through centerline, in units of 104Ur/Lr). (b) Top and (c) side views of isosurfaces of

vorticity magnitude (|ω| = 3.1×10−7 Ur/Lr). Motion of virtual wake-creating body is from

left to right. Views shown represent full streamwise, but subregions of lateral and vertical,

domain (table 1).

15
3.2 Mean defect history

The history of the maximum mean velocity defect Ud is shown in figure 7a and b. The Ud

variation is smoother than that typically found from simulations, due to the reduced statistical oscillations associated with the larger-than-usual streamwise domain used here (the

same can be said for the mean-wake height and width histories in figure 7c, d). In contrast

to the unstratified-wake DNS of RCC, for which the difference between the maximum and

the centroid of the mean defect increased with time, and becomes significant at later stages,

here there is essentially no difference between the two, for the time period considered. There

is however a noticeable difference in the location in the y-z plane at which the two occur;

as in RCC, we use the location of the centroid of the mean-velocity distribution u(y, z) to

define the current/local wake centerline coordinates (yc, zc).

Figure 7b also shows histories of the streamwise u

, lateral v and vertical w velocity

fluctuations at the wake centerline (symbols). These, along with the Ud history, reveal that

it is not until after N t ≈ 1 that the wake turbulence is fully developed. We therefore view

the early phase, from N t = 0 until N t = N tvo ≈ 1, as effectively a precursor simulation. The

exact virtual-time origin tvo is defined by the unstratified, 3D similarity diagnostic U −3/2

d

shown in figure 7a, from which we deduce N tvo = 0.645 (chain-double-dotted line). In order

to remove some of the uncertainty associated with comparing the present results during the

3D and NEQ regimes with those from experiments and other simulations, the histories that

follow will be presented as functions of Nt

 ̃= N t−N tvo = N t−0.645. Figure 7a also indicates

that Ud continues to exhibit pure-3D behavior until Nt

 ̃ = N t − N tvo ≈ 4.8 − 0.645 ≈ 4.2,

since this is when the similarity diagnostic first deviates from the linear variation.

The early wake is well represented by the canonical Ud ∼ t−2/3 history despite not being

fully self similar during this time: although fully turbulent and essentially unstratified and

axisymmetric (figure 8a), the wake’s non-dimensional Reynolds-shear-stress profile in cylindrical coordinates −v

xvr/U2

d (not shown) is not solely a function of the non-dimensional

radius r/h (i.e. with no time dependence). This is consistent with the evolution of the unstratified, vortex-ring-initialized, Case VR wake described in RCC, for which (non-universal)

self-similarity did not appear until after the time corresponding to N t ≈ 11 for the present

simulation – and universal self-similarity not until Nt > 70. (For the present DNS, the

relationship between N t and the non-dimensional time used in RCC is t/t∗ = 30N t.) In

fact, in the range corresponding to 1 < Nt < 11, the Ud history for Case VR decays as t

−1,

because of the unique turbulence distribution induced by the Case VR rings (which produce

an eddy viscosity profile that is temporally constant and radially uniform; see figures 8a and

9b of RCC). Here the ring initialization, which used somewhat larger streamwise spacing,

2.1ho rather than Case VR’s 1.85ho, is such that Ud ∼ t−2/3 is a better representation than

t

−1 for the 3D regime. These observations illustrate the profound impact the initial wake

structure can have upon its development.

As observed elsewhere (e.g. in Spedding’s [1] experiments and Gourlay et al.’s [13]

simulation), the magnitude of the rate of decay |dUd/dt| slows markedly for Nt > 5, as

the flow enters the NEQ phase, such that Ud ∼ t−n with n ≤ 1/4. In this region, the Ud

variation is well represented by the interpolant 0.25(N t + 1)−1/4. (It should be mentioned

that power laws with slightly different exponents could also reproduce the Ud decay to good

accuracy, given the flexibility afforded by the use of a finite virtual origin for time.) The

time at which the 3D and NEQ power laws intersect, Nt

 ̃ = N t − N tvo ≈ 6 − 0.645 ≈ 5.5,

can be viewed as an estimate of the boundary between the two regimes for the present flow.

The generality of this finding will be discussed in § 3.6.

The Ud history for the variable-viscosity run (Case R1500; see § 2) begins to significantly

differ from the constant-viscosity/variable-Rd result after Nt

 ̃ ≈ 100 (compare dotted and

solid lines in inset plot in figure 7b), which implies non-inertial effects are important after

this time. The NEQ period revealed by the evolution of Ud begins at Nt

 ̃ ≈ .5.5 and

16
     	 
  

	



	



N t

(Ud/Udo)−3/2

(a)

     







 







N t

Ud/Udo

−0.75

−0.85

(b)

    



	







  





Nt

 ̃= N t − N tvo

z/ho

1/2

βz

(c)

    



	



  

	



	







Nt

 ̃

y/ho

1/4

0.35

βy

(d)

Figure 7. Histories of (a, b) maximum mean velocity Ud, (c) mean height z and vertical

spreading rate βz and (d) mean width y and lateral spreading rate βy: , present

DNS (with z = r∗

z and y = r∗y in c and d, respectively);  (in b only), 0.1(u

u)

1/2 c /Ud

(arbitrary factor of 0.1 is to aid visibility); Δ (in b only), 0.1(vv)

1/2 c /Ud;  (in b only),

0.1(w

w)

1/2 c /Ud; (in c and d), hz and hy; (in c and d), r

z and ry; ◦ (in d only),

L

σ/r∗y (cf. Meunier & Spedding [17]); • (in d only), Lσ/r∗y (cf. Meunier & Spedding [17]); ,

(in c) βz = (1/Ud)dr∗

z /dt and (in d) βy = (1/Ud)dr∗y/dt; , (in c) βz = (1/Ud)dr∗z /dt and (in

d) βy = (1/Ud)dr∗

y/dt for Case R1500 (variable-viscocity (Ud r∗/ν ≥ 1500) run); ,

curve fits for (a, b) Ud, (c) r∗

z and (d) r∗y, where in (a) (Ud/Ud(t))−3/2 = 3.1(N t−0.645), in (b)

Ud/Udo = 0.47(N t−0.645)−2/3 for Nt < 10, Ud/Udo = 0.25(N t+1.0)−1/4 for 3 < Nt < 100,

Ud/Udo = 2.05(N t−10)−0.85 for Nt > 40, while in (c) and (d) r∗

z /ho = r∗y/ho = 1.23( Nt

 ̃)1/3

for Nt < 10 where Nt

 ̃= N t−N tvo and N tvo = 0.645 (see figure 21b). ⊕ (in b only), time at

which Case R1500 first reaches its minimum Reynolds number, Ud r∗/ν = 1500; (inset

plots in b, c, d only), Case R1500 (with z = r∗

z and y = r∗y in c and d). At N t = 0,

r∗

z = r∗y = r∗o = 0.865ho. The Ud(t) interpolants intersect at N t ≈ 6 and 45.

17
follows (or remains just above) the NEQ power-law curve fit until about Nt

 ̃ = 30 and

Nt

 ̃= 50 respectively for the constant- and variable-viscosity (i.e. lower and higher Reynolds

number) cases. The effect of viscosity on the NEQ duration is even more apparent in

the height histories (figure 7c). The Reynolds numbers for both cases are thus too low to

allow an extended NEQ region – and that the period during which significant turbulence

is ‘regenerated’ by Kelvin-Helmholtz instabilities (Diamessis et al. 2011), if it exists at all,

will be quite brief (cf. figure 5b). (see instantaneous lateral vorticity contours at Nt

 ̃= 20 in

figure 5 of Redford et al. 2014). The Reynolds-number dependence of the duration of the

NEQ regime will be considered further in § 3.6. For now we note that the Ud history seen

here during the final stage of the NEQ phase is more consistent with the t

−1/2 behavior of the

‘buoyancy-controlled’ subregion between NEQ and Q2D proposed by Meunier et al. (2006)

than the Ud plateau of the ‘NEQ2’ subregion identified by Brucker & Sarkar (2010).

A t

−3/4 power law for Ud could be assumed for the Q2D regime (cf. [3,25]), which would

agree well with the present data for N t ≥ 80 if a virtual origin of N t = 30 were used.

(Recall the earlier comment about the interplay between the power law’s exponent and

virtual origin.) However, a better fit is provided by t

−0.85 for N t ≥ 60, with a smaller virtual

origin (N t = 10) (chain-double-dotted curve in figure 7a). The −0.85 exponent is also more

consistent with conservation of mean momentum, when coupled with the behavior of the

vertical and lateral defect dimensions (see below). The NEQ/Q2D boundary defined by the

intersection of the corresponding power-law interpolants is Nt

 ̃≈ 55; the t

−0.85 interpolant

closely corresponds to the Ud variation until the end of the (primary) simulation, N t ≈ 840.

Note that the final time is not large enough for a purely viscous 3D regime to appear (such

that Ud ∼ t−1; see Meunier et al. [3]). In fact, the TKE remains a significant fraction of

U2

d for all the times considered here (figure 15a inset plot). Nevertheless, the importance of

viscous effects during the Q2D phase is underscored by the significant deviation for Nt > 90

of the Ud histories from the constant- and variable-viscosity flows (figure 7b inset plot).

The height histories are shown in figure 7c. Both the ‘half-radius’ height hz (i.e. half the

distance between the locations at which the mean velocity in the vertical plane through the

centerline, u(z,yc), is 0.5 of Ud; not a best-fit to an assumed Gaussian profile), and an integral measure r∗

z are shown; the latter is defined as r∗z = [(2 ln 2/Ud)

∞

0 u(rz, yc)rzdrz]

1/2,

where r2

z = (z − zc)2 and (yc, zc) are the coordinates of the local wake centerline. The

integral height has been normalized such that it is equal to hz when u(z,yc) is Gaussian.

Consequently, the greater the difference between hz (chain-dotted line) and r∗

z (solid) in

figure 7c, the more non-Gaussian the vertical mean velocity profile (figure 9a).

During the 3D phase, when the wake is effectively unstratified, its height and width

grow in proportion to t

 ̃1/3, and its non-dimensional mean growth/spreading rate βz =

(1/Ud)dr∗

z /dt is approximately constant. This can be seen in figure 7c, for Nt < ̃ 3 (chaindotted curve and symbols, which respectively trace a t

 ̃1/3 interpolant for r∗

z and the βz

variation). Perhaps unsurprisingly, given the flattening and widening effect the stratification

has on the mean defect, the height/width interpolant deviates from the DNS sooner (near

N t − N tvo = 2) than the corresponding Ud power-law interpolant does (Nt

 ̃≈ 4.5).

Despite the wake not becoming fully self-similar, the spreading rate during the unstratified regime remains fixed near 0.03, which is roughly half that observed at much later

times (corresponding to Nt > 70) for the self-similar growth of RCC’s unstratified Case VR

wake. The streamwise velocity fluctuations at the centerline are also weaker here during

the 3D regime, than for Case VR during its (first, non-universal) fully self-similar phase,

with ( u

u )

1/2 c /Ud ≈ 0.2 here (open squares in figure 7b) compared to approximately 0.3 for

Case VR.

The difference between the integral height for the constant- and variable-viscosity runs

(solid and dotted curves, respectively, shown in inset plot), and that of the non-dimensional

mean growth/spreading rate βz (solid and open diamonds) – namely that the region of zero

18
slope (and thus zero βz) lasts longer for larger Rd (Nt

 ̃ ≈ 70 rather than 30) – corroborate

Meunier et al.’s [3] prediction, and Diamessis et al.’s [2] observation, that the duration of

the NEQ regime increases with increasing Reynolds number.

Figure 7c also shows the differences in the evolution of the height of the mean defect

and vertical extent of the wake turbulence, which become apparent when one compares

histories of hz (chain-dotted curve) or r∗

z (solid) with that of rz (dashed), a TKE based

integral length, defined as r

z = 


{k (z − zc)2}/{k}

1/2 where k = 1

2u

iui. The mean defect

does indeed ‘collapse’ early in the NEQ regime, in the sense that both hz and r∗

z decrease

with time from Nt

 ̃≈ 8 until about Nt ̃= 15 – although the reduction in the integral height

is much weaker than that exhibited by the u = 0.5Ud measure hz. As found by Gourlay et

al. [13], the turbulent region continues to spread vertically during the NEQ regime, when

dr∗/dt ≈ 0, before an even more dramatic collapse occurs, with dr

z/dt < 0 (much larger

than the negative dhz/dt and dr∗

z /dt seen earlier) between Nt

 ̃≈ 20 and 100. This behavior

was also observed for the self-propelled case by Lin & Pao [10] in their flow-visualisation

experiments, and in Brucker & Sarkar’s [9] DNS of self-propelled and towed wakes.

After the wake enters the Q2D regime, r

z begins to grow again, as do hz and r∗z . All

three heights then converge, as the vertical spreading becomes a purely viscous phenomenon,

such that hz ≈ r∗

z ≈ rz [13] and each grow as t

1/2 [3]. This can be seen in figure 8, which

includes contours of the mean defect (lines) and TKE (colour) at times representing the

three regimes. The tendency for the height of the turbulent region to occupy a smaller and

smaller fraction of the mean defect as time passes is clear. (The changes to the TKE budget

corresponding to the change in shape of the TKE profiles – for example, the growth of the

‘flattened, double-pronged’ structure in the y-z plane that begins to appear during the NEQ

phase – will be considered in § 3.5.) The non-dimensional vertical growth rate βz is therefore

positive but much smaller than during the 3D regime, when the spreading is effected by the

turbulence.

Due to the conservation of mean streamwise momentum, the product of the mean height,

width and maximum velocity will remain approximately constant through the wake’s lifetime. (It would be exactly constant for a pure 2D-Gaussian mean defect. The Udhzhy

product for the primary simulation remains within ±16% of its value at N t = 1 for the

full 1 ≤ N t ≤ 840 range.) The variations of the power-law exponents for Ud and hz during

the 3D, NEQ and Q2D regimes are thus expected to lead to variations for the mean-width

history (which would conflict with Diamessis et al.’s [2] claim that hy ∼ t1/3 for all three

regimes). For example, the combination of monotonically decreasing Ud and dr∗

z /dt → 0

requires the mean width to grow faster than it does during the 3D phase as the wake enters the NEQ regime. Figure 7d illustrates the enhanced lateral growth, compared to the

3D power law (chain-double-dotted curve), in terms of both the mean half-width/radius

hy (chained-dotted curve) and the integral width r∗

y = [(2 ln 2/Ud)

∞

0 u(zc, ry)rydry]

1/2,

with r2

y = (y − yc)2 (solid curve), for 1.5 ≤ Nt

 ̃ ≤ 4.5. Thereafter, the rate of spreading

slows, falling from a maximum near Nt

 ̃ = 5 to a local near-zero minimum near Nt ̃ = 20

(see open and closed diamonds in figure 7d, which show the history of the non-dimensional

lateral spreading rate βy = (1/Ud)dr∗

y/dt). Beginning from Nt

 ̃ ≈ 25, as the wakes leaves

the NEQ and enters the Q2D regime (and Ud(t) begins to increase its rate of decay), the

lateral growth increases rapidly until it reaches a value comparable to that observed during

the unstratified 3D period, with r∗

y ∼ t0.35. This is similar to the hy ∼ t1/3 proposed by

Diamessis et al. [2], but differs from the t

1/4 included in Meunier et al.’s [3] model, and the

t

0.23 observed in Brucker & Sarkar’s [9] DNS, for the Q2D range. The difference between

the n = 0.35 and 0.25 exponents is shown by the thin-solid lines in figure 7d, which supports

our preference for the larger value – which in turn explains (in light of the viscous n = 1/2

growth of hz, and conservation of momentum) our use of n = −0.85 to fit the Ud(t) history

(figure 7b). (Meunier et al.’s [3] exponent combination of (−3/4, 1/2, 1/4) for (Ud, hz, hy)

19
      
 












 






(y − y0)/ho

(z − z0)/ho

(a) Nt

 ̃= 1.2

0.06

0.03

0       
 












 

	



(y − y0)/ho

(b) Nt

 ̃= 2.5

0.07

0.035

0

    



	



	

 






(z − z0)/ho

(c) Nt

 ̃= 18

      










 

	



(z − z0)/ho

(d) Nt

 ̃= 82

 	  	 



	



	

 	

	



(y − y0)/ho

(z − z0)/ho

(e) Nt

 ̃= 586

Figure 8. Mean-defect and TKE contours at: (a) Nt

 ̃ = 1.2, (b) Nt ̃ = 2.5, (c) Nt ̃ = 18,

(d) Nt

 ̃= 82 and (e) Nt ̃= 586. TKE (shaded contours) normalized by local U2

d . Solid-line

contours correspond to u = 0.05Ud and 0.5Ud.

20
also satisfies the constant-momentum constraint.) Despite the profound qualitative differences in the turbulence structure in the 3D and Q2D regimes, the rate at which the mean

defect spreads laterally is about the same.

The width of the turbulence, like that of the mean defect, grows continuously throughout

the lifetime of the wake. This can be seen by the history of the TKE-based integral width

r

y = 


{k (y − yc)2}/{k}

1/2 in figure 7d (dashed line). In contrast to the height of the

turbulent region, its width remains significantly larger than the width of the mean defect,

with r

y/hy > 1 for each of the 3D, NEQ and Q2D periods (cf. Gourlay et al. [13]). (We

cannot say whether or not the waviness in r

y(t) for Nt > ̃ 20 is physical.) The closed and open

circles in figure 7d respectively correspond to the ratios L

σ/r∗y and Lσ/r∗y, where Lσ and Lσ

are, respectively, the distance from the wake centerline to the maximum values of the lateral

Reynolds shear stress −u

v and the horizontal-plane velocity fluctuations (uu + vv)1/2.

Integral version of these scales are used, based on the profiles presented in Meunier &

Spedding [17], with L

σ = ∞

−∞ |u

v|dy/u

2

0 and Lσ = ∞

−∞(u

u + vv)1/2dy/(2π)1/2u



0 ,

where u



2

0 = exp(0.5)(u

v)max, u



0 = 0.5 exp(0.5)(u

u + vv)

1/2 max, and (u

v)max and (uu +

vv)

1/2 max are the maxima over y of |u

v| and (uu + vv)1/2, respectively.

In contrast to r

y/r∗y, which is always significantly greater than one, the shear-stress width

ratio L

σ/r∗y remain close to unity for the entire run, slightly less before Nt

 ̃ ≈ 10, slightly

more thereafter, while for the horizontal TKE measure 0.7 < L

σ/r∗y ≤ 1 throughout. The

L

σ ≈ r∗y and Lσ ≈ r∗y behavior, which is similar to that exhibited by the experimental data

for various bluff bodies given in Meunier & Spedding [17], implies a type of self-similarity

for the lateral profiles, especially for −u

v, that survives as the importance of stratification

and viscous effects grows, and the flow moves through the 3D, NEQ and Q2D regimes.

The decoupling of the vertical and lateral spreading processes, and the associated greater

importance for the latter of the mixing due to turbulence relative to molecular viscosity,

is illustrated by the much closer agreement of the mean width r∗

y from the constant- and

variable-viscosity flows (figure 7d inset) compared to the corresponding r∗

z histories (figure 7c

inset).

Figure 9 contains vertical (9a) and lateral (9b) profiles of streamwise-mean velocity

u, normalized by the local maximum Ud and height hz or width hy of the defect. The

better collapse and smoother profiles seen here, compared to those from Gourlay et al.’s

[13] DNS and Diamessis et al.’s [2] LES, are thought to be a symptom of the better

averaging sample afforded by the present larger streamwise domain. Included in these plots

are results from times correspond to the structures shown in figures 4 (dashed curve; 3D

regime), 5 (chain-dashed; NEQ) and 6 (long-dashed; Q2D). (The dotted profile is the mean

induced by the vortex ring initialization.) Given the qualitative differences in the 3D, NEQ

and Q2D structures, the collapse at these times is striking, as is the agreement with the

Gaussian idealization (open symbols), particularly near the centerline. At the earlier times,

the non-Gaussian behavior away from the centerline may be an artifact of the vortex-ring

initialization, and/or due to attenuation of the eddy viscosity, associated with corrugations

of the interface between the vortical and irrotational regions [31] (figure 4).

For the lateral profiles, which tend to be more consistently Gaussian than the vertical

ones, the near-centerline agreement is similar to that found by Spedding et al. [25]. (The

unphysical lateral asymmetry at Nt

 ̃ = 137 is evidence of finite eddy sample at this time.

Curiously, the u(zc, y) profile at Nt

 ̃ = 586, which is the result of averaging over an even

smaller number of eddy structures than are present at Nt

 ̃= 137, is nearly exactly symmetric;

see long-dashed curve in figure 9b.) The vertical profiles, on the other hand, can at times

exhibit significant non-Gaussian features, across the entire wake. One might suspect the

‘extra peakiness’ near the centerline during the early stages of the Q2D regime (Nt

 ̃ = 82

and 137) is due to the coarsening of the vertical resolution across the mean defect (i.e.

Δz/hz = 0.14; cf. table 1) mentioned at the end of § 2. However, this possibility was ruled

21
      



	



(z − zc)/hz

u(y0, z)/Ud (a)

      



	



(y − yc)/hy

u(y, z0)/Ud (b)

Figure 9. Profile of mean streamwise velocity in (a) vertical and (b) lateral planes through

the wake centerline: , N t = 0; , Nt

 ̃ = 1.2 , Nt ̃ = 2.5; , Nt ̃ =

20; , Nt

 ̃ = 82; , Nt ̃ = 137 , Nt ̃ = 586; ◦, Gaussian idealisation

u/Ud = exp(− ln(2) s2), where s is either (y − yc)/hy or (z − zc)/hz. Coordinates (y0, z0)

and (yc, zc) respectively define the initial and current location of the wake centerline, given

by the centroid of the mean velocity defect.

out by repeating this portion of the simulation with increased Mz, such that Δz/hz = 0.07,

resulting in u(z,yc) variations that are nearly indistinguishable to those shown in figure 9a.

The mechanisms responsible for the evolving shape of the mean defect revealed in figure 9,

as well as its magnitude, height and width (figure 7), can be further understood by examining

the corresponding evolution of the terms in the mean-momentum transport equation. But

before attending to that task, we consider the history and budget of mean kinetic energy.

3.3 Integrated mean kinetic energy history and budget

The effect of the stably stratified background upon the three components of mean kinetic

energy (MKE) is illustrated in figure 10a. Histories of the (cross-domain) area-integrated

quantities {Kαα} = { 1

2uαuα} for α = 1, 2, 3 (no sum on α) are shown. Only during the

NEQ regime are the lateral and vertical energies non-negligible compared to the streamwise

component, and even then {K11} remains at least 1.5 orders of magnitude larger than {K22}

and {K33}. After Nt

 ̃≈ 3, the buoyancy suppresses the vertical component compared to the

lateral one, and there is a net transfer of energy from the vertical mean kinetic energy K33

to the mean potential energy (MPE), 1

2 g{ρρ}/ρ∞|Γ|.

This can be seen in figure 10b, which documents the variation of the terms in the integrated mean-kinetic-energy transport equation:

d{Kαα}

dt = −{Pααjj}−{Dααjj} + {Bαα}, α = (1, 2, 3), (5)

where

−Pααjj = u

αu

j

∂uα

∂xj

is the mean-to-turbulence transfer term (i.e. the negative rate of TKE production),

Dααjj = ν

∂uα

∂xj

∂uα

∂xj

22
    









Nt

 ̃

{Kαα}/{Kii}o

(a)

    





	







Nt

 ̃

−

d{K11}/dt budget

(b)

Figure 10. (a) Components of integrated mean kinetic energy {Kαα} = { 1

2uαuα} (no sum

on α): , {K11}; , {K22}; , {K33}. Values normalized by MKE at N t = 0,

{Kii}o = {K11}o+{K22}o+{K33}o. (b) Terms in the budget of the streamwise mean kinetic

energy {K11}: , −{D1122} = −{ν(∂u/∂y)2}; , −{D1133} = −{ν(∂u/∂z)2};

, −{D1122}−{D1133}; , −{P1122} = {uv ∂u/∂y}; , −{P1133} =

{u

w ∂u/∂z}; , −{P1122}−{P1133}; , −{D1122}−{D1133}−{P1122}−{P1133};

, {B33} = −g{ρ w}/ρ∞. Data for −g{ρ w}/ρ∞ for Nt > ̃ 8 and −{P1133} for Nt > ̃ 18

(when both quantities are less than 10−7) not included due to excessive statistical scatter.

Budget terms given in units of U3

r Lr.

is the mean dissipation rate, and

Bαα = −guαρ δα3/ρ∞

is the mean buoyancy transfer (which appears with opposite sign in the transport equation for the mean potential energy 1

2 gρ ρ/ρ∞|Γ|). Note that the integral of both the

advection Aααjj = −uj∂Kα/∂xj and flux-divergence/transfer Tααjj = −∂uαu

αu

j/∂xj +

ν ∂[uα∂uα/∂xj ]/∂xj terms are zero, {Aααjj} = {Tααjj} = 0. It will be useful to distinguish

between mean stratification effects associated with the vertical and lateral gradients. For the

present flow, for which ∂( )/∂x1 = ∂( )/∂x = 0, we can write −Pααjj = −Pαα22 − Pαα33

and Dααjj = Dαα22 + Dαα33, where −Pαα22 = u

αv ∂uα/∂y, −Pαα33 = uαw ∂uα/∂z,

Dαα22 = ν(∂uα/∂y)2 and Dαα33 = ν(∂uα/∂z)2.

The open symbols in figure 10b trace the evolution of the buoyancy term. Its magnitude

and sign (i.e. {B33} = −g{w ρ}/ρ∞ < 0) imply that the buoyancy term acts as a pure

sink of MKE, and source of MPE, which is largest early in the NEQ regime, until at least

Nt

 ̃≈ 9 (it is not shown at later times because of its statistical uncertainty, exemplified by

rapid small-amplitude oscillations of both signs). There are thus no mean ‘re-stratification

effects’ associated with re-conversion of mean-potential to mean-kinetic energy (cf. [2]) for

the present flow’s NEQ phase. The possibility that potential energy supplies kinetic energy

to the turbulence during the NEQ regime will be re-examined below.

The end of the 3D period is first evident in the mean dissipation {D1122} and {D1133}, in

that the mean vertical gradients enhanced by the vertical flattening of the wake overwhelm

the lateral variations, thereby arresting the reduction of {D1133} (compare long- and shortdashed curves in figure 10b). Because the dissipation terms are much smaller than the

MKE-to-TKE transfer terms, their effect on the net d{K11}/dt is not as pronounced as

that caused by the precipitous fall in magnitude of the vertical component {−P1133}, and

the recovery (and eventual reversal in slope) of the lateral component {−P1122}, that begin

23
ï

ï

ï

ï











(z − z0)/ho

(a) ×10−4

−10 0 10

−10

0

10

∂u/∂t

ï

ï

ï

ï











(b) ×10−4

−10 0 10

−10

0

10

−v ∂u/∂y − w ∂u/∂z

ï

ï

ï

ï











(z − z0)/ho

(y − y0)/ho

(c) ×10−4

−10 0 10

−10

0

10

−∂uv/∂y − ∂uw/∂z

ï

ï

ï

ï











(y − y0)/ho

(d) ×10−4

−10 0 10

−10

0

10

ν (∂2u/∂y2 + ∂2u/∂z2)

Figure 11. Terms in the mean streamwise momentum budget at Nt

 ̃= 2.5: (a) rate of change

∂u/∂t; (b) advection −v ∂u/∂y − w ∂u/∂z; (c) Reynolds-stress divergence −∂u

v/∂y −

∂u

w/∂z; (d) viscous-stress divergence ν (∂2u/∂y2 + ∂2u/∂z2). Terms given in units of

U2

r /Lr, where Ur = 2.09Udo and Lr = ho.

near Nt

 ̃= 11. Brucker & Sarkar [9] found that the rapid attenuation of P1133 is due to the

stratification-induced collapse of the −u

w Reynolds stress. The terms combine to yield a

near-plateau of d{K11}/dt between Nt

 ̃= 20 and 50, before the flow enters its Q2D phase,

when the net MKE decay is roughly split between the vertical dissipation and lateral TKE

production terms.

3.4 RANS budget

In this section we examine the terms in the budget of mean streamwise momentum, and

how stratification modifies their spatial structure in vertical (y-z) planes normal to the wake

axis, to better understand the mechanisms responsible for the decay and spread of the mean

defect during the 3D, NEQ and Q2D regimes. For the present time-dependent parallel flow

the streamwise Reynolds-averaged Navier-Stokes (RANS) equation can be written

∂u

∂t =



−v

∂ u

∂y − w

∂ u

∂z 

  

advection

+



−∂ u

v

∂y − ∂ u

w

∂z 

  

Reynolds-stress divergence

+ ν

∂2 u

∂y2 +

∂2 u

∂z2



  

viscous-stress divergence

. (6)

During the 3D phase, the reduction in Ud and increase in the wake height and width

24
– i.e. the negative ∂u/∂t near the centerline and the positive ∂u/∂t at the wake edge –

is driven by the divergence of the Reynolds stress, −∂u

v/∂y − ∂uw/∂z. This can be

inferred from figure 11, which shows the cross-wake plane distribution of ∂u/∂t and how

the advection, turbulence and viscous terms contribute to it at Nt

 ̃= 2.5. At this time, the

axisymmetric/3D nature of the flow has been slightly altered by the buoyancy, as the flow

moves towards the NEQ regime. Although the relative strengths of the vertical and lateral

components of the Reynolds-stress divergence remain comparable at the edges of the wake,

the flattening and widening of the region of wake turbulence caused by the buoyant forcing

is apparent (figure 11c). The reduced rate of vertical spreading is due to vertical advection,

with −w ∂u/∂z < 0 above and below the centerline (dark regions in figure 11b), while the

lateral spreading is enhanced by the corresponding advection term, in that −v ∂u/∂y > 0

at the left and right of the defect (light regions in figure 11b).

Within the NEQ regime the story is more complex. Figure 12a shows that at Nt

 ̃ =

18, when dhz/dt ≈ dr∗

z /dt ≈ 0, the regions at the top and bottom of the wake contain

negligible ∂u/∂t, while u continues to grow at the sides. Significant internal wave activity,

in the κx = 0 (i.e. streamwise mean) mode, is generated by the wake, revealed by the

vertical mean momentum flux divergence, −∂u w/∂z (figure 12b). However, the waves do

not make a meaningful contribution to the net momentum (or the kinetic energy), since

−∂u w/∂z − ∂u v/∂y ≡ −w∂u/∂z − v∂u/∂y is nearly identically zero outside the wake

(figures 12c and d). The (slower, compared to 3D) rate of decrease of Ud(t) is now driven

by lateral Reynolds-stress divergence −∂u

v/∂y and vertical advection −w∂u/∂z, while

the lateral component −v∂u/∂z retards the lateral spreading near the centerline plane but

increases it further above and below (figures 12d). The Reynolds stresses, especially the

vertical component −u

w, are profoundly affected by the buoyancy (which drives w → 0;

see figure 15b), such that the importance of their divergence has been reduced to levels

comparable to that of the vertical viscous stresses. (Contrast figures 12e, f and g). That

viscous effects are important within the NEQ regime was illustrated by the difference in the

defect histories for the primary simulation and Case R1500 (figure 7). The lateral advection

and the vertical viscous terms (which tend to mirror and counteract each other) reveal a

‘flattened-X’ pattern in the mean velocity that is not apparent in the vertical u(z,yc) and

lateral u(zc, y) profiles in figure 9. The lateral viscous-stress divergence (not shown) remains

negligible.

The difference between the ∂u w/∂z contours in figures 12b (Nt

 ̃= 18) and 13b (Nt ̃= 82)

displays the tendency for the strength of the internal waves to fade as the flow passes from

the NEQ to Q2D regime. Once the wake enters the Q2D phase, the core deceleration is

dominated by the lateral Reynolds-stress divergence, and to a lesser extent by the vertical advection and vertical viscous-stress terms. This can be seen in figures 13e, c and f,

respectively. The vertical Reynolds-stress −∂u

w/∂z and lateral viscous-stress ν ∂2u/∂y2

terms are now negligible. We also note that the multi-layered structure first seen in the

instantaneous vorticity contours earlier in the NEQ range, at Nt

 ̃ ≈ 20 (figure 5), has now

matured to the point it affects the RANS budget, leaving its signature in the terms involving spatial gradients of the mean velocity, such as −v∂u/∂y (figure 13d). At the top and

bottom edges of the wake, the vertical viscous-stress divergence is now solely responsible

for the small net positive ∂u/∂t here – a fact foreshadowed by the t

1/2 growth of the mean

height observed in figure 7c. In contrast, the (faster) lateral spreading is due to the lateral

advection (figure 13d) and especially the lateral Reynolds-stress term (figure 13e).

As time passes, the importance of the −∂u

v/∂y and ν ∂2u/∂z2 terms continues to

grow relative to the other ones, until they are the only significant contributors to the RANS

budget (cf. Meunier et al. [3]). At Nt

 ̃ = 586 (figure 14), both decelerate the core region

by about the same amount, but have much different effects on the spreading of the defect:

the mean acceleration at the top and bottom of the defect caused by the viscous diffusion is

much weaker than that produced at the sides by the lateral divergence of the −u

v Reynolds

25
(z − z0)/ho

(a)

∂u/∂t

−20 0 20

−10

0

10

×10−5

1.5

1

0.5

0

−0.5

−1

−1.5

(z − z0)/ho

(b)

−∂u w/∂z

−20 0 20

−10

0

10

×10−5

1.5

1

0.5

0

−0.5

−1

−1.5

(z − z0)/ho

(c)

−w ∂u/∂z

−20 0 20

−10

0

10

×10−5

1.5

1

0.5

0

−0.5

−1

−1.5

(z − z0)/ho

(d)

−v ∂u/∂y

−20 0 20

−10

0

10

×10−5

1.5

1

0.5

0

−0.5

−1

−1.5

(z − z0)/ho

(e)

−∂uw/∂z

−20 0 20

−10

0

10

×10−5

1.5

1

0.5

0

−0.5

−1

−1.5

(z − z0)/ho

(f)

−∂uv/∂y

−20 0 20

−10

0

10

×10−5

1.5

1

0.5

0

−0.5

−1

−1.5

(z − z0)/ho

(y − y0)/ho

(g)

ν ∂2u/∂z2

−20 0 20

−10

0

10

×10−5

1.5

1

0.5

0

−0.5

−1

−1.5

Figure 12. Terms in mean streamwise momentum budget at Nt

 ̃= 18 (NEQ regime): (a) rate

of change ∂u/∂t; (b) −∂u w/∂z; (c) vertical advection −w ∂u/∂z; (d) lateral advection

−v ∂u/∂y; (e) vertical Reynolds-stress divergence −∂u

w/∂z; (f) lateral Reynolds-stress

divergence −∂u

v/∂y; (g) vertical viscous-stress divergence ν ∂2u/∂z2. Lateral viscousstress divergence ν ∂2u/∂y2 is negligible compared to terms shown. Terms given in units of

U2

r /Lr, where Ur = 2.09Udo and Lr = ho.

26
(z − z0)/ho

(a)

∂u/∂t

−60 0 60 −20

0

20

×10−6

3

2

1

0

−1

−2

−3

(z − z0)/ho

(b)

−∂u w/∂z

−60 0 60

−20

0

20 ×10−6

3

2

1

0

−1

−2

−3

(z − z0)/ho

(c)

−w ∂u/∂z

−60 0 60 −20

0

20

×10−6

3

2

1

0

−1

−2

−3

(z − z0)/ho

(d)

−v ∂u/∂y

−60 0 60

−20

0

20 ×10−6

3

2

1

0

−1

−2

−3

(z − z0)/ho

(e)

−∂uv/∂y

−60 0 60 −20

0

20

×10−6

3

2

1

0

−1

−2

−3

(z − z0)/ho

(y − y0)/ho

(f)

ν ∂2u/∂z2

−60 0 60

−20

0

20 ×10−6

3

2

1

0

−1

−2

−3

Figure 13. Terms in mean streamwise momentum budget at Nt

 ̃= 82 (Q2D regime): (a) rate

of change ∂u/∂t; (b) −∂u w/∂z; (c) vertical advection −w ∂u/∂z; (d) lateral advection

−v ∂u/∂y; (e) lateral Reynolds-stress divergence −∂u

v/∂y; (f) vertical viscous-stress divergence ν ∂2u/∂z2. Vertical Reynolds-stress divergence −∂u

w/∂z and lateral viscousstress divergence ν ∂2u/∂y2 are negligible compared to terms shown. Terms given in units

of U2

r /Lr, where Ur = 2.09Udo and Lr = ho.

27
(z − z0)/ho

(a)

∂u/∂t

−100 0 100

−30

0

30

×10−8

5

0

−5

(z − z0)/ho

(b)

−∂uv/∂y

−100 0 100

−30

0

30

×10−8

5

0

−5

ï

(z − z0)/ho

(y − y0)/ho

(c)

ν ∂2u/∂z2

−100 0 100

−30

0

30

×10−8

5

0

−5

Figure 14. Terms in mean streamwise momentum budget at Nt

 ̃ = 586 (late Q2D regime):

(a) rate of change ∂u/∂t; (b) lateral Reynolds-stress divergence −∂u

v/∂y; (c) vertical

viscous-stress divergence ν ∂2u/∂z2. Vertical mean-flow momentum divergence −∂u w/∂z,

advection −v ∂u/∂y − w ∂u/∂z, vertical Reynolds-stress divergence −∂u

w/∂z and lateral

viscous-stress divergence ν ∂2u/∂y2 are negligible compared to terms shown. Terms given

in units of U2

r /Lr, where Ur = 2.09Udo and Lr = ho.

stress, induced by the pancake-eddy turbulence (cf. figure 6). (The lateral asymmetry of the

−∂u

v/∂y contours in figure 14b is a statistical anomaly, due to the relatively small number

of eddy structures included in the streamwise average at this time.) The decoupling of the

viscous/vertical and turbulent/lateral mixing mechanisms during the late-Q2D regime is

consistent with the disparate rates of vertical and lateral spreading observed in figures 7(c, d),

with βz larger than βy for Nt > ̃ 100, and the earlier observation that the growth of hz is a

purely viscous phenomenon during this time.

3.5 Turbulence kinetic energy history and budget

Figure 15a summarizes the manner in which stratification alters the cross-wake-averaged

potential and kinetic energy of the turbulence, 1

2 gρ

ρ/ρ∞|Γ| (symbols) and k = 1

2 u

iui

(lines), respectively. (Recall the cross-wake average is over regions where u ≥ 0.05 Ud.)

The density fluctuations are largest early in the NEQ regime, when buoyancy most

strongly affects the flow, but contain negligible energy compared to that in the velocity

fluctuations throughout the wake’s lifetime (note the factor of 103 multiplying the potential

energy in figure 15a).

The present flow, with its much larger initial Froude number (F∞ ≈ 130), does not contain the buoyancy-induced temporal oscillations of TKE observed at early times in Dommermuth et al.’s [8] F∞ = 4 LES and Brucker & Sarkar’s [9] F∞ = 8 DNS, for which dk/dt

first becomes positive near Nt

 ̃= 5.5 for the former and Nt ̃= 2.3 for the latter (see figure 21

regarding values of the virtual origin Nt

 ̃vo used for those simulations). Those oscillations,

which are damped more slowly in the lower F∞ flow, are indicative of transfer back and

forth between the potential and kinetic energy of the turbulence, associated with changes in

sign of the TKE-to-turbulence potential energy (TPE) exchange term, Btke = −gρ

w/ρ∞

28
   





	









    

	









Nt

 ̃

103 × 1

2 ρρ/(Ud|Γ|/N)2, k/U2do

k/U2

d

−4/3

(a)

   



	





	



	



Nt

 ̃

3bij 3bij

(b)

Figure 15. (a) Cross-wake-averaged turbulence potential and kinetic energy histories: ◦, 1

2 ρρ

; , k; , 1

2 u

u; , 1

2 vv

; , 1

2 ww. Inset plot shows

components scaled by local Ud. (b) Reynolds stress anistropy bij = |u

iuj|/2k − δij/3

histories: , b11; , b22; , b33; , b12; , b13; , b23.

(Brucker & Sarkar mention, but do not quantify, ‘large values of the buoyancy flux at early

time’). For the present flow, Btke always acts as a net sink of TKE (figure 3b) despite the

monotonic decrease with time of the local mean Froude number Fd (figure 1a, b).

Instead of increasing towards the beginning of the NEQ regime, here the TKE increases

for 30 < Nt < ̃ 70, i.e. during the NEQ-to-Q2D transition defined by the Ud history (figure 7b). The reasons for this increase, which are unrelated to the TPE-to-TKE transfer via

Btke, will soon be clear.

The difference between the rates of decay of the three velocity fluctuations, as well as the

non-constant variation of k/U2

d for Nt > ̃ 4 (solid curve in figure 15a inset), demonstrate

the anisotropy introduced by the stratification, and how it breaks the early-time axisymmetry, leading to a flow for which the relationship between the magnitude of the velocity

fluctuations and the maximum mean defect continues to evolve throughout the NEQ and

Q2D regimes. Of the three components, the streamwise TKE 1

2 u

u (dashed curve) exhibits the least variation with respect to the local U2

d , but even so varies significantly (by

over a factor of five) as the wake regimes come and go. For Nt > ̃ 150, when the TKE is

reasonably approximated by a t

−4/3 power law, k decays more slowly than U2

d does, such

that (assuming Ud ∼ t

−0.85 in the Q2D range) k/U2

d ∼ t+0.37 (figure 15a inset plot).

The disparate effect of the stratification upon the three velocity-fluctuation components,

and thus on the structure of the Reynolds-stress tensor, can also be seen in figure 15b.

These Reynolds-stress-anisotropy histories show that the velocity fluctuations first deviate

from axisymmetric-decay behavior (i.e. vv ≈ ww and u

v ≈ uw) near Nt

 ̃= 2, the time

at which the mean-defect height and width begin to differ from each other (figure 7c, d).

Thereafter, the vertical component falls monotonically, first slowly then precipitously, such

that ww and both u

w and vw are negligible for Nt > ̃ 40 (cf. Spedding [12] and Diamessis

et al. [2]). The absence of vertical velocity means that buoyancy can no longer directly cause

a transfer between the turbulence potential and kinetic energies.

The streamwise normal stress uu becomes a larger and larger fraction of the 2k =

uu + vv + ww total until the flow is well within the NEQ regime, when dk/dt becomes

positive, near Nt

 ̃ ≈ 30. Afterwards, the normal component vv increases faster than u

u,

and u

v is the only significant shear stress. The dominance of uv corresponds to that of

the P1122 = −u

v∂u/∂y production term; recall the importance of −{P1122} in the {K11}

29
     








	







Nt

 ̃

d{k}/dt budget

−1.7

(a)







    





Nt

 ̃

{Ptke}/{tke} 3{ ̃αα}/{tke}

(b)

Figure 16. Terms in the budget of integrated TKE {k}: , rate of production {Ptke};

, rate of dissipation {tke}; , (minus) rate of buoyancy transfer −{Btke}; (in

b only), {P1122}; thin-solid ( ) (in b only), {P1133}; (in b only), {P2222};

(in b only), {P3333};  (in b only), { ̃11}; ◦ (in b only), { ̃22}; • (in b only), { ̃33};

(in a only), A(Nt

 ̃)−5/3 = A(N t − N tvo)−5/3, with (A,Ntvo) = (0.00138, 0.645).

Production terms not shown are either zero because of homogeneity ({P1111} = {P2211} =

{P3311} ≡ 0) or negligible ({P2233}≈{P3322} ≈ 0). Components of dissipation rate

 ̃αα = ν(∂u

α/∂xj)(∂uα/∂xj ) (no sum on α), where  ̃ii = tke − ν∂2u

iuj/∂xi∂xj , and

tke = 2νs

ij sij . Terms in (a) given in units of U3

r Lr, in (b) as fraction of local {tke}.

budget and see below. During the Nt > ̃ 70 phase of the Q2D regime, when dk/dt is again

negative (and the flow is dominated by pancake-eddy structures), the lateral fluctuations

have more energy than the streamwise fluctuations do, with v

v/uu → 2 for Nt > ̃ 400.

We now consider the transport equation for the TKE k,

∂k

∂t = Atke + Ptke + Ftke + Btke − tke + Wtke, (7)

where the rates of change due to advection Atke, production Ptke, flux divergence Ftke,

buoyancy transport Btke and dissipation tke are, respectively, Atke = −uj∂k/∂xj, Ptke =

Piijj = −u

iuj∂ui/∂xj , Ftke = ∂fj/∂xj , with fj = −1

2u

juiui − ujp

/ρ∞ + ν∂u

iuj/∂xi +

ν∂k/∂xj, and Btke = −gρ

u

iδi3/ρ∞ and tke = 2νsij sij , with s

ij = 1

2 (∂u

i/∂xj + ∂uj/∂xi).

For the flow at hand, these reduce to Atke = −v∂k/∂y − w∂k/∂z, Ptke = Pii22 + Pii33,

Ftke = ∂[−1

2 vu

iui−v

p/ρ∞+ν∂vv/∂y+ν∂vw/∂z+ν∂k/∂y]/∂y+∂[−1

2wu

iui−wp

/ρ∞+

ν∂v

w/∂y +ν∂ww/∂z +ν∂k/∂z]/∂z and Btke = −gρw/ρ∞. The work done by the fringe

forcing Fi in (1) is Wtke = u

iFi (i = 1, 2, 3).

We begin by examining the cross-plane-area-integrated histories of the terms in (7),

shown in figure 16. (This information was presented in a different format in figure 3b;

it is repeated here to highlight the relative importance of individual budget terms.) The

cross-plane averages used in figure 15 can be approximated from the cross-plane integrations

used in figure 16 by dividing the latter by the cross-sectional area of the mean defect, such

that e.g. k≈{k}/4πhyhz. Since the advection and flux-divergence terms integrate to

zero, figure 16 reveals the contributions of the pure source Ptke, sink tke and TKE-TPE

transfer Btke. The cross-plane variation of these terms, and Atke(y, z) and Ftke(y, z), will be

investigated below.

For the present high-F∞ flow, the fringe-work contribution {Wtke} to the TKE balance

is negligible (recall that the production, dissipation and buoyancy sum shown in figure 3b is

30
   











 





(z − zc)/hz

(a) Nt

 ̃= 2.5

   











 







(b) Nt

 ̃= 5.5

   











 







(z − zc)/hz

(c) Nt

 ̃= 18

   











 	



	

(d) Nt

 ̃= 40

   











 







(y − yc)/hy

(z − zc)/hz

(e) Nt

 ̃= 82

   











 





(y − yc)/hy

(f) Nt

 ̃= 586

Figure 17. Contours of buoyancy transfer Btke as a fraction of (tke)max, the local maximum

tke from figure 20, at: (a) Nt

 ̃ = 2.5, (b) Nt ̃ = 5.5, (c) Nt ̃ = 18, (d) Nt ̃ = 40, (e)

Nt

 ̃ = 82 and (f) Nt ̃ = 586. Minimum (darkest) contour level -Clim corresponds to the

most negative Btke/(tke)max at each time. White solid-line contours indicate five levels of

positive Btke/(tke)max, over range +0.01|Clim|≤Btke/(tke)max ≤ [Btke/(tke)max]max, where

[Btke/(tke)max]max is maximum Btke/(tke)max, which at Nt

 ̃ = 2.5, 5.5, 18, 40, 82 and 586

is respectively +0.0007, +0.085, +0.37, +0.14, +0.06 and +0.005. Black solid-line contours

correspond to u = 0.05Ud and 0.5Ud.

31
   











 





(z − zc)/hz

(a) Nt

 ̃= 2.5

   











 





(b) Nt

 ̃= 5.5

   











 







(z − zc)/hz

(c) Nt

 ̃= 18

   











 	


	




(d) Nt

 ̃= 40

   











 








(y − yc)/hy

(z − zc)/hz

(e) Nt

 ̃= 82

   











 





(y − yc)/hy

(f) Nt

 ̃= 586

Figure 18. Contours of density-fluctuation variance ρρ at: (a) Nt

 ̃ = 2.5, (b) Nt ̃ = 5.5,

(c) Nt

 ̃ = 18, (d) Nt ̃ = 40, (e) Nt ̃ = 82 and (f) Nt ̃ = 586. Values given in units of

(ρ∞U2

r /gLr)2. Solid-line contours correspond to u = 0.05Ud and 0.5Ud.

32
sufficient to provide an excellent representation of the d{k}/dt history, for 0 ≤ Nt < 840.)

The net buoyancy transfer is also much smaller than the production and dissipation, but is

not negligible, with −{Btke} reaching a maximum of about 30% of the integrated dissipation

{tke} from early in the NEQ regime (Nt

 ̃ ≈ 5) and remaining fixed near this ratio until

Nt

 ̃ ≈ 70 (figure 16b), after which −{Btke} is a monotonically decreasing fraction of {tke}

until it becomes negligible after Nt

 ̃≈ 150. Between Nt ̃≈ 25 and 70, all three non-negligible

terms are roughly constant in absolute terms, with d{tke}/dt and d{Btke}/dt approximately

zero and d{Ptke}/dt increasing slightly. Spedding [12] found a similar NEQ plauteau in the

total dissipation tke + Diijj in his high-F∞ experiments, which he attributed to ‘vertical

shear taking over’. We have already observed that the global/integrated buoyancy transfer

{Btke} acts as a pure sink throughout the lifetime of the present flow (figure 3b). At times,

however, there are local regions in which Btke is positive, where the potential energy of

the turbulence is converted to kinetic energy. During the early NEQ phase, these tend

to be concentrated just above and below the wake (see white-line contours in figure 17b),

while later they are more diffuse and disorganized, although still occurring mostly outside

the wake core (figure 17c-e). Much later, in the viscous-dominated Q2D regime, when

{Btke}/{tke} → 0, positive Btke tends to lie in the non-vortical regions to the left and right

sides of the wake (figure 17f). The relationship between Btke and the density fluctuations

ρρ can be inferred by comparing figures 17 and 18. The regions of maximum negative Btke

are in general highly correlated with those of maximum |ρ|, although the departure from

axisymmetry associated with the beginning of the NEQ regime is more apparent in the ρ

ρ

field at Nt

 ̃= 2.5 (figure 18a) than in the corresponding Btke distribution (figure 17a). The

local ρ

ρ maxima at Nt

 ̃ = 2.5 at the top and bottom of the wake, near z = ±1.5hz, are

consistent with the ‘excess’ TPE that later ‘springs back’ into TKE at these locations, by

way of the positive Btke observed in figure 17b at Nt

 ̃= 5.5.

Figure 16a shows that the Reynolds number is large enough for the early/3D-regime

evolution to satisfy the inertial scaling tke ∼ k3/2/h, such that tke ∼ t−7/3 and {tke} ∼

t

−5/3 (chain-double-dot curve; cf. Brucker & Sarkar [9]). For late times, Nt > ̃ 200 or

so, the decay of {Ptke} and especially {tke} is reasonably well approximated by t

−n with

n = 1.7 ≈ 5/3. Thus, paradoxically, the TKE dissipation evolves according to the inertial

scaling during both the 3D and Q2D regimes, despite their qualitatively different eddy

structures (see figure20d).

Once the wake turbulence is fully established, the net rate of TKE dissipation is larger

than that of the production Ptke, except for a substantial period during the NEQ regime,

from Nt

 ̃≈ 20 to 90, when {Ptke} is significantly larger than {tke} (and gradually increases,

in absolute terms), after which it returns to its pre-NEQ behavior, with {Ptke}/{tke} < 1 for

the Q2D phase. This explains the growth of TKE during this period, mentioned above. From

figure 16b we find that the enhanced NEQ-regime production is due solely to amplification

of the P1122 = −u

v ∂u/∂y component (dotted curve), since the −uv stress grows more

rapidly than ∂u/∂y decreases (and the fall of u

w is more rapid than the attenuation

of the decay of ∂u/∂z; see figure 16b and histories of the MKE dissipation components

D1122 = ν(∂u/∂y)2 and D1133 = ν(∂u/∂z)2 in figure 10). This ‘structural’ alteration,

represented by the decrease of |u

w| and increase of |uv|, was also found by Brucker &

Sarkar [9] in their DNS for both towed and self-propelled cases. They concluded that

stratification causes the TKE production to mainly consist of P1122 during this time. The

present results underline the extent to which this is true, and that in fact Ptke can effectively

be replaced by P1122 for all times after Nt

 ̃≈ 20.

While both the present and Brucker & Sarkar simulations contain enhanced production,

only here does it lead to growth of TKE (and not just reduction in the net TKE decay)

during the NEQ regime. This is because of the large amount of energy radiated by internal

waves between Nt

 ̃= 25 and 75 in Brucker & Sarkar’s F∞ = 8 flow, which counteracts the

impact of the P1122 growth. The strength of the internal waves is an indication of when

33
the TKE growth occurs. When they are strong, at low F∞, buoyancy leads to an early

TKE increase driven by the TPE-to-TKE transfer, before the NEQ regime begins – during

which the energy the waves radiate mitigates the effect of the enhanced TKE production

associated with the flattened, two-layer NEQ structure, and d{k}/dt is small but negative.

When the waves are weak, at high F∞, they do not attenuate the enhanced Ptke ≈ P1122,

leading to net TKE growth during the NEQ regime.

Also included in figure 16b are histories of the three components of the integrated rate of

homogeneous dissipation, { ̃ii} = { ̃11}+{ ̃22}+{ ̃33}, as a fraction of the full (actual, thermodynamic) rate of TKE dissipation {tke} = 2ν{s

ijsij}, where  ̃ii = ν(∂ui/∂xj )(∂ui/∂xj ).

(There is very little difference between the integrated homogeneous and full/inhomogeneous

dissipation for this flow, with 0.9999965 ≤ { ̃ii}/{tke} ≤ 1.0012 for all times considered.)

Examining  ̃αα for α = 1, 2, 3 allows us to assess the interplay between the rates of production and dissipation of the individual normal-stress components u

αuα. (Other dissipation

decompositions have been utilized by Spedding [12] and Brucker & Sarkar [9]; see also

Itsweire et al. [32].)

The vertical component { ̃33}∼{(∂w

/∂x)2} + {(∂w/∂y)2} + {(∂w/∂z)2} falls monotonically as w does, with d{ ̃33}/dt becoming negative from Nt

 ̃ ≈ 7 until Nt > ̃ 50, when

{ ̃

33} is negligible with respect to {tke}.

For the streamwise component, the increase in P1122 for Nt > ̃ 5 is accompanied by an

increase in { ̃11}, such that the net effect (to which the pressure-strain correlation φ11 =

p(∂u/∂x)/ρ∞ may also contribute) is to arrest the decay of uu, such that d{uu}/dt first

changes sign and then becomes slightly positive as the flow enters its NEQ phase (figure 15a).

This is an example of what is often referred to as an ‘extra strain’ effect (Bradshaw [33]),

which involves a difficult-to-model change that results from the imbalance of much larger

terms of opposite sign, each of which change much more than the net imbalance. Similar

behavior is found in the Q2D regime, in that as d{u

u}/dt becomes negative again, both

{P1122} and { ̃11} decrease (the latter accounting for roughly one-third of the total {tke}).

The evolution of the lateral component during the 3D-to-NEQ transition is not directly

affected by production (since P2211 ≡ 0 and P2222 ≈ P2233 ≈ 0). Consequently, the slight

increase in {v

v} between Nt

 ̃ ≈ 25 and 80 shown in figure 15a must be due to intercomponent transfer via the pressure-strain term φ22 = p(∂v

/∂y)/ρ∞, mitigated by { ̃22},

which first falls slightly below {tke}/3 (i.e. the isotropic value approximated by all three

components before Nt

 ̃≈ 5), then grows to account for about 2/3 of the total.

The dissipation histories are broadly consistent with those observed in Brucker & Sarkar’s

[9] low-F∞/self-propelled wake simulation, in that (∂u

/∂z)2 and (∂v/∂z)2 (i.e. the dominant terms in  ̃11 and  ̃22, respectively) are respectively the second-largest and largest terms

for Nt > ̃ 10 (although the contributions of ∂v

/∂x, ∂w/∂x and ∂w/∂z were not quantified, these are likely to be negligible). The { ̃αα} histories are thus also consistent with

the enhanced vertical gradients associated with the flattened eddy structures of the NEQ

and Q2D regimes (figures 5 and 6). We note that { ̃22}/{ ̃11} ≈ 2 for Nt > ̃ 150 while

v

v/uu > 1.7 in the late-Q2D regime, i.e. for Nt > ̃ 500 (figure 15b).

Terms in the TKE budget for the (essentially) unstratified axisymmetric benchmark at

Nt

 ̃ = 1.2 are illustrated in figure 19. Shown are radial profiles of the production Ptke,

dissipation tke, buoyancy transport Btke (which at this time is very small) and the flux

divergence Ftke, which solely redistributes TKE within (and thus integrates to zero over)

the y-z plane, and thus played no role in the integrated budget histories considered to

this point. The importance of the flux divergence (open symbols) is indicated by the large

contribution it makes to the net ∂k/∂t (solid symbols), especially to the TKE growth at

the edge of the wake. The shapes of the production and dissipation profiles are the same as

those found in the self-similar axisymmetric experiment of Uberoi & Freymouth [34], but

here the transport increases the TKE not just at the edge of the wake (as it does in the

experiment) but also near the centerline. (The difference is perhaps a consequence of all the

34
      







±r/h

∂k/∂t budget

Figure 19. Radial profiles of terms in the TKE budget at Nt

 ̃= 1.2: , rate of production

Ptke; , (minus) rate of dissipation −tke; , rate of buoyancy transfer Btke; ◦, flux

divergence Ftke; , Ptke − tke + Btke + Ftke. Advection −v∂k/∂y − w∂k/∂z is negligible

compared to terms shown. Budget terms normalized by U3

r /Lr, where reference velocity

Ur = 2.09Udo and reference length Lr = ho.

fluxes, including the viscous ones, being included in Ftke in figure 19, rather than only the

inviscid, boundary-layer terms shown in [34].) But while the profile shapes of the budget

terms for the Uberio & Freymuth experiment and the present simulation are comparable,

their relative magnitudes are not – which provides further evidence for the non-universal

nature of the turbulence for (unstratified) wakes generated by different objects [5,35],RCC.

The influence of the stratification upon the three dominant terms in the TKE budget,

and in particular how they evolve in the y-z plane, can be observed in figure 20. To facilitate

comparison, ‘quadrant-averaged’ contours of ∂k/∂t, Ptke, −k and Ftke at selected times are

combined in each plot, in the ‘northwest’, ‘northeast’, ‘southwest’ and ‘southeast’ quadrants,

respectively. The lateral and vertical coordinates are scaled by their respective width or

height, to clarify the spatial variations, and reveal the early ‘kinematic’ stratification effects,

which act solely to flatten and widen the axisymmetric distributions. The black lines trace

the shape and extent of the mean defect, via 0.05Ud and 0.5Ud contours, while the white-line

contours indicated the dominant term(s) at each time.

At Nt

 ̃= 2.5 (figure 20a), the only indication of the beginning of the end of the 3D regime

is the slight lack of axisymmetry in the flux term, which causes the net TKE growth at the

edge of the wake to be larger at the sides than at the top and bottom. The production

and dissipation are not yet deeply affected, since using the width and height to normalise

the cross-plane axes causes them to maintain their (axisymmetric) radial shapes from the

unstratified state (figure 19). The location of the maximum production is still just inside the

u = 0.5Ud (i.e. r = h) contour (close to the location of the maximum TKE; see figure 8b),

and the maximum dissipation falls between r = h and the centerline. As is the case for

the unstratified/axisymmetric state, the flux term plays a large role in setting the net

∂k/∂t variation, and is solely responsible for the TKE growth at the wake edge (since the

dissipation tends to cancel the production here).

35
   











 	

	



(z − zc)/hz

(a) Nt

 ̃= 2.5

∂k/∂t Ptke

−tke Ftke

   











 





(b) Nt

 ̃= 18

∂k/∂t Ptke

−tke Ftke

   











 	

	



(y − yc)/hy

(z − zc)/hz

(c) Nt

 ̃= 82

∂k/∂t Ptke

−tke Ftke

   











 







(y − yc)/hy

(d) Nt

 ̃= 586

∂k/∂t Ptke

−tke Ftke

Figure 20. Contours of terms in TKE budget at: (a) Nt

 ̃= 2.5, (b) Nt ̃= 18, (c) Nt ̃= 82 and

(d) Nt

 ̃= 586. Streamwise-averaged quantities have been further averaged (‘double folded’)

across the lateral y and vertical z axes through the wake centerline (yc, zc), and plotted in

separate quadrants, such that the upper left (‘northwest’), upper right (‘northeast’), lower

left (‘southwest’) and lower right (‘southeast’) region of each plot respectively corresponds

to the net rate of change ∂k/∂t = Atke + Ptke + Ftke + Btke − tke, production Ptke, (minus)

dissipation −tke and flux-gradient transport Ftke. Buoyancy flux term Btke is shown in

figure 17. Advection Atke = −v∂k/∂y − w∂k/∂z is negligible compared to the other terms.

Budget terms normalized by U3

r /Lr, where reference velocity Ur = 2.09Udo and length Lr =

ho. White solid-line contours indicate ±0.67 of the larger of the minimum and maximum

magnitude at each time. Black solid-line contours correspond to u = 0.05Ud and 0.5Ud.

The (maximum,minimum) values at Nt

 ̃ = 2.5, 18, 82 and 586 are respectively (+1.49 ×

10−5, −2.25 × 10−5), (+1.77 × 10−7, −1.1 × 10−7), (+5.3 × 10−8, −3.2 × 10−8) and (+8.15 ×

10−11, −2.6 × 10−10).

36
When the flow is well into the NEQ regime, Ftke is much less significant, leaving the

production-dissipation imbalance to control ∂k/∂t (figure 20b, for Nt

 ̃= 18). The dissipation

has by this time become quite anisotropic, with tke ≈ 2 ̃11 ≈ 4 ̃22 ≈ 5 ̃33 (figure 16b).

Nevertheless, its cross-plane variation – while tending to align into two flat regions above

and below the y axis (consistent with the double-layer structure observed in figure 5, for

Nt

 ̃ ≈ 20) – remains more axisymmetric than that of the production. In contrast, Ptke

tends to be concentrated into a ‘wing-like’ region near y = ±hy, because it is now nearly

all associated with the −u

v stress (Meunier et al. [3]), with Ptke ≈ P1122 = −uv∂u/∂y.

Consequently, the net ∂k/∂t ≈ Ptke − tke tends to be segregated according to sign, into

either an ‘hourglass’-shaped zone aligned with the z axis (TKE decay) or double ‘hotspots’

of TKE growth focused near y = ±hy. The latter is compatible with the ‘double-pronged’

structure in the TKE contours at Nt

 ̃≈ 20, seen in figure 8c.

The flow’s entrance into its Q2D phase is characterized by the continuing importance of

Ptke ≈ P1122, which is again concentrated along the y axis, near y = ±hy (figure 20c, for

Nt

 ̃ = 82). However, the production now contains a small region of negative Ptke centered

at the origin, due to negative P2222 = −v

v∂v/∂y. This is due to the combination of

finite v

v (whose maximum occurs at the centerline), and the positive ∂v/∂y associated

with the mean stagnation flow (i.e. vertical collapse and lateral expansion) at the origin.

The dissipation is now even more concentrated in two flat regions above and below the

y axis, within the u = 0.5Ud contour, than it was at Nt

 ̃ = 18 – corresponding to the

enhanced vertical shear mentioned by e.g. Spedding [12], and the primary roles played at

this time by  ̃11 ∼ (∂u

/∂z)2 and  ̃22 ∼ (∂v/∂z)2 (figure 16b). The flux-divergence term

has grown in importance, compared to its Nt

 ̃= 18 levels, such that it significantly reduces

Ptke (near y = ±hy) and decreases tke (near z = ±hz). The net result leads to the TKE

distribution shown in figure 8d, which demonstrates the persistence of the flattened doublepeaked structure as the flow passes from the NEQ to the Q2D regime (see also Brucker &

Sarkar [9]).

Much later, at times for which the turbulence is dominated by the quasi-2D pancake

eddies displayed in figure 6, the Ptke ≈ P1122 production is still active near y = ±hy, and tke

(≈ 1.5 ̃22) is still focused in horizontal layers near z = ±0.5hz, but both are weaker than the

flux term Ftke. As is the case for Ptke and tke, the distribution of Ftke is closely correlated

with the two-layer pancake-eddy structure at this time. The largest (positive) verticalgradient contribution ∂f3/∂z is aligned with the horizontal layers centered near z = ±hz/2

(and thus reduces the effect of the dissipation here), while the smallest (negative) values lie

between the layers, along the y axis; the maximum and minimum lateral divergence ∂f2/∂y

are clustered along the y axis, at the lateral edges of the dual Ptke (and TKE) peaks, with

the negative values closer to the wake centerline. The flux term combines with Ptke and tke

to yield large negative ∂k/∂t in a shallow horizontal layer near the centerline. While there

is still some TKE growth due to the combined effect of Ptke and Ftke, between the y and z

axes along the u = Ud contour, it is very weak.

3.6 Comparison with previous simulations

Figure 21 shows the mean-defect histories from earlier simulations at various Froude and

Reynolds numbers, along with the present results (Fdo = 130, Rdo = 2625). Available

DNS and a recent LES are included: The DNS of Gourlay et al. [13] (Fdo = 34, Rdo =

1472) and Brucker & Sarkar’s [9] Case TR50F04 (Fdo ≈ 0.75, Rdo ≈ 3250) and TR10F20

(Fdo ≈ 3.75, Rdo ≈ 650), and Diamessis et al.’s [2] Case R100F64 LES (Fdo ≈ 34, nominal

Rdo ≈ 11, 000). Also shown are the unstratified/3D-regime power-law diagnostics for each

case (inset plots in figure 21a, b). These are used to estimate virtual-time origins for each

simulation, to facilitate comparison by plotting histories in terms of Nt

 ̃= N t − N tvo.

The various histories cannot be expected to agree in the scaling used for figure 21, which

37
     









 	 





 	 

	



 	 

	



 	 





 	 







N t

N t

N t

N t

N t

Nt

 ̃

Ud/Udo

χ

χ

χ

(a)

    







 	 





 	 

	



	

 	 

	



	

 	 

	

 	 

	



	

N t

Nt Nt

N t

Nt

 ̃

z/y

ξ

ξ ξ

(b)

    





Nt

 ̃

z/o

(c)

    





Nt

 ̃

y/o

(d)

Figure 21. Histories of (a) maximum velocity Ud, (b) height-to-width ratio z/y, (c) height

z and (d) width y of mean defect: , present; ◦, [13]; , Case TR50F04 of [9]; •,

Case TR10F20 of [9]; , Case R100F64 LES of [2] . Inset plots in (a) and (b) show unstratified

similarity diagnostics for (a) velocity χ = (Ud/Udo)−3/2 and (b) height ξz = (z/o)3 (lower

curve/symbols) and width ξy = (y/o)3 (upper curve/symbols), which are used to define the

virtual-time origin tvo for each case, with N tvo = 0.645 for the present DNS, N tvo = 0.2 for

Brucker & Sarkar’s TR50F04, N tvo = −0.5 for Diamessis et al. and N tvo = 0 for Gourlay et

al. [13] and Brucker & Sarkar’s [9] TR10F20. For the present DNS, the linear curve-fit for χ

and ξ ( in upper right-hand-side inset plots in a and b) are χ = (Ud/Udo)−3/2 = B · Nt

 ̃

and ξ = (r∗/r∗

o)3 = C · Nt

 ̃, where Nt ̃ = N t − N tvo and (B, C, N tvo) = (3.1, 2.9, 0.645).

The maximum defect velocity Ud in (a) for Brucker & Sarkar’s [9] Case TR50F04 (but not

TR10F20) is replaced by their Um.k.e.

1 . Integral measures are used for the height and width

for the present DNS, with z = r∗

z , y = r∗y, and o the value of r∗z = r∗y at N t = 0. For

Brucker & Sarkar’s [9] Case TR50F04, z and y are respectively represented by integral

scales (in their notation) R3 and R2, with o the half-radius of the initial Gaussain meanvelocity profile, while for Diamessis et al. [2], and effectively for Gourlay et al. [13], z = hz,

y = hy and o = ho, where hz and hy are determined from fitting an assumed Gaussian

profile.

38
takes no account of Froude and/or Reynolds numbers, since the maximum velocities Ud,

heights z, and widths y are presented solely in terms of their initial values. (Also, z is

for some cases given in terms of the 0.5Ud half-height hz and for others in terms of integral

measures; see figure 21c caption.) Nevertheless, several common features can be observed,

including similar power-law behavior over the same Nt

 ̃ range. The height decreases with

time for Nt > ̃ 2 in every case, although the reduction of hz is larger than that experienced

by the integral-measure heights (which is consistent with the behavior found for the present

DNS, in that both dhz/dt and dr∗

z /dt are negative for 8 < Nt < ̃ 18, but hz falls faster than

r∗

z ; figure 7b).

The 3D-to-NEQ transition also tends to occur at about the same time for each of the

simulations, with Nt

 ̃ between 2 to 2.5 the point at which the the heights and widths first

begin to differ from each other; this can be seen in figure 21b, in both the height-to-width

ratio z/y and is even more obvious in the unstratified/axisymmetric height/width diagnostic ξ = (/o)3 shown in the inset plots (recall that linear variation with time of ξ is

characteristic of the pure-axisymmetric 3D regime).

The reasonably good collapse of the mean aspect ratio z/y over the early NEQ periods

of the five simulations was not anticipated, both because of the wide Froude and Reynolds

number ranges, as well as the differences in the initialization strategies used for the various

simulations (e.g. Brucker & Sarkar’s [9] DNS held the mean velocity fixed until N t =

0.2 (TR10F20) and 1 (TR50F04), while Diamessis et al.’s [2] LES gradually introduced

the background stratification during a precursor simulation). The z/y histories again

demonstrate the increase of the Nt

 ̃ range of the NEQ regime with Reynolds number –

compare the times at which the z/y histories begin to deviate from the negative powerlaw behavior for NEQ, for Brucker & Sarkar’s [9] TR10F20 (40 < Nt < ̃ 100), Gourlay et

al. [13] (Nt

 ̃≈ 70), present (Nt ̃≈ 105), and Brucker & Sarkar’s [9] TR50F04 (Nt > ̃ 400, at

least). (The late-time z/y history for Diamessis et al. [2] (open diamonds) is thought to be

statistically unreliable, in light of the relatively small streamwise domain used.) The same

message regarding the Reynolds-number dependence of the NEQ period can be inferred

from figure 21c, by examining the durations of the dz/dt ≈ 0 regions for each case.

The tendency for z and y to deviate from the 3D/axisymmetric ideal sooner than

Ud does, which was seen above (figure 7), also holds for the other simulations (although

not as markedly), since the unstratified/axisymmetric similarity diagnostic for Ud, χ =

(Ud/Udo)−3/2, remains linear longer (until Nt

 ̃≈ 2 to 2.5) than the corresponding height/width

diagnostic, ξ (Nt

 ̃≈ 2 to 5); cf. figures 21a and b inset plots.

A scaling to account for the significant Froude-number sensitivity of the Ud histories

seen in figure 21a has been developed by Meunier & Spedding [17]. It amounts to assuming

Ud will eventually depend solely on the volume-flux deficit Id, the Brunt-V ̈ais ̈al ̈a frequency

N and time t, which leads via dimensional analysis to

Ud

(IdN2)1/3 = Φ(N t − N t0), (8)

where Φ is a universal non-dimensional function of N t − N t0 and N t0 is a non-universal

offset time (unrelated to the 3D virtual-origin time N tvo), which will depend for example

on Reynolds number and the near-field/early-time wake evolution. In terms of a virtual

wake-creating sphere of diameter D emersed in a steady uniform stream of density ρ∞

and speed V∞, experiencing a drag force FD, we can write Id = 1

8V∞D2CD, where CD =

FD/(πρ∞V 2

∞D2/8). Equation (8) then implies

Ud

V∞

F2/3 ∞ = f(CD)Φ(N t − N t0), (9)

where f(CD)=(πCD/2)1/3. Equation (9) shows the limitation of attempting to collapse

Ud/V∞ data from wakes generated by various bluff bodies (with variable CD) by using only

39
 	  	  	 



	



	



	



N t

(Ud/I1/3

d N2/3)−4/3

(a)

   







N t − N t0

Ud/I1/3

d N2/3

−0.75

(b)

 	  	  	 



	



	



N t

(Ud/I1/3

d N2/3)−1/0.85

(a)

   







N t − N t0

Ud/I1/3

d N2/3

−0.85

(b)

Figure 22. Histories of rescaled mean defect (a, c) as Q2D power-law/similarity diagnostic

and (b, d) with respect to Q2D time shift N t0, assuming (a, b) n = 3/4 and (c, d) n = 0.85:

, present DNS; ◦, Gourlay et al. [13]; , Case TR50F04 of Brucker & Sarkar [9];

•, Case TR10F20 of [9]; , Case R100F64 LES of Diamessis et al. [2]. Thin-solid curves

( ) in (a, c) are linear curve fits of Q2D diagnostic (which for [13] and Case TR10F20

of [9] are nearly coincident), used to define time shift N t0 in (b, d), such that for (b) (n = 3/4)

N t0 =: 19 (present DNS), −55 ( [13]), −63.5 (Case TR50F04 of [9]), −57 (Case TR10F20

of [9]) & 11 (Case R100F64 of [2]), and for (b) (n = 3/4) N t0 =: 7 (present DNS), −70

( [13]), −66.5 (Case TR50F04 of [9]), −71 (Case TR10F20 of [9]) & −3 (Case R100F64

of [2]).

40
F∞ (e.g. by plotting (Ud/V∞)F2/3 ∞ versus N t; cf. [1,9]). Following Meunier & Spedding [17],

(8) can be recast in terms of the virtual effective diameter Deff = D(CD/2)1/2 for bluff bodies

of arbitrary shape,

Ud

V∞

F2/3

eff = π1/3Φ(N t − N t0), (10)

where Feff = 2V∞/NDeff is the effective Froude number.

The validity of (8) is tested in figure 22. Meunier & Spedding [17] found the universality

under consideration cannot be expected until the flow reaches the Q2D regime (see their

figure 5b). Consequently, we assume Φ(N t − N t0) = CQ2D(N t − N t0)−n, where CQ2D is

a universal constant and n ≈ 3/4. Similar slopes of the (Ud/I1/3

d N2/3)−1/n diagnostic are

found for all five simulations, for both n = 3/4 (figure 22a, with CQ2D = 0.08−3/4 ≈ 6.65)

and the n = 0.85 suggested by the present DNS (figure 22c; CQ2D = 0.055−0.85 ≈ 11.8).

Note the dependence of the beginning of the region of common slope (straight thin-solid

lines in figure 22a, c), and the offset times N t0 (N t-axis intercept), on the Reynolds number

and initialization scheme used by each simulation. The reasonable collapse of the results in

figures 22b and d, produced by using the N t0 defined for each case by the corresponding

curve fits in figures 22a and c, points to the validity of (8), and thus to that of the founding

assumptions used in the Meunier & Spedding [17] analysis.

4 Summary and closing comments

The DNS results presented above reveal the effect of a weakly stratified, buoyantly stable

background upon the structure and statistics of a turbulent wake, during its various phases

of development, including its 3D, NEQ and Q2D regimes. This has been documented by

the evolution of the mean momentum and energy, turbulence kinetic energy, and the terms

in their budgets. The dependence of the quality of the statistics on adequate ‘eddy sample

size’ has been observed, as has the Reynolds-number dependence of the duration, in terms

of N t, of the NEQ regime.

The mean momentum budget clarifies the relationship between the turbulence suppression due to buoyancy and the turbulent and viscous mechanisms responsible for the vertical

and lateral growth of the wake. The TKE budget, on the other hand, quantifies the manner in which the changes in turbulence structure are reflected in changes to the relative

importance and composition of the terms whose imbalance defines the net rate of change

of TKE. Surprisingly, despite these changes in structure, the TKE dissipation decays at

approximately the same rate during the early/3D and late/Q2D regimes. This affirms yet

again the dissimilarity between the Q2D state and pure 2D turbulence [11].

The present results illustrate that the wake collapse, signaling the beginning of the

NEQ regime, need not be accompanied by the transfer from the potential energy of the

turbulence to the kinetic energy of the turbulence, as has been proposed elsewhere. Instead,

the enhanced TKE exhibited during the NEQ regime is solely the result of buoyancy-induced

structural changes, which lead to increased TKE production. Thus, while the potential-tokinetic energy transfer is observed at early times for lower F∞ cases (Dommermuth et al. [8];

Brucker & Sarkar [9]), it is not an essential feature of the 3D to NEQ transition.

For the present high-F∞ case, internal wave activity is most pronounced at the start

of the NEQ regime. However, even during this time, internal waves make a negligible

contribution to the turbulence energy budget. This is in contrast to the strongly stratified

cases studied elsewhere, for which the amount of turbulence energy radiated by internal

waves during the NEQ regime can be comparable to that produced by the mean shear or

lost via viscous dissipation (cf. Brucker & Sarkar [9]).

The role of stable stratification in driving the wake towards a universal state has also

been displayed, in that the Q2D structure is independent of the details of the wake initializa41
tion (a fact frequently observed in other experiments and simulations). The present study

also indicates that during the Q2D period, the history of the maximum mean defect Ud (appropriately normalized, in terms of the initial volume-flux deficit Id and the Brunt-V ̈ais ̈al ̈a

frequency N) is a unique function of nondimensional time N t (relative to a flow-dependent

virtual offset/origin). This supports the validity of the framework proposed by Meunier &

Spedding [17] in their universal model of the mean-defect evolution for a stably stratified

towed wake.

References

1. Spedding, G. R.: The evolution of initially turbulent bluff-body wakes at high internal

Froude number. J. Fluid Mech., vol. 337, 1997, pp. 283–301.

2. Diamessis, P. J.; Spedding, G. R.; and Domaradzki, J. A.: Similarity scaling and vorticity structure in high-Reynolds-number stably stratified turbulent wakes. J. Fluid Mech.,

vol. 671, 2011, pp. 52–95.

3. Meunier, P.; Diamessis, P. J.; and Spedding, G. R.: Self-preservation in stratified momentum wakes. Phys. Fluids, vol. 18, 2006, p. 106601.

4. Tennekes, H.; and Lumley, J.: A First Course in Turbulence. The MIT Press, 1972.

5. Bevilaqua, P. M.; and Lykoudis, P. S.: Turbulence memory in self-preserving wakes. J.

Fluid Mech., vol. 89, no. 3, 1978, pp. 589–606.

6. Chomaz, J. M.; Bonneton, P.; Butet, A.; and Hopfinger, E. J.: Vertical diffusion of

the far wake of a sphere moving in a stratified fluid. Phys. Fluids, vol. 5, no. 11, 1993,

pp. 2799–2806.

7. Riley, J. J.; and Lelong, M.-P.: Fluid motions in the presence of strong stratification.

Annu. Rev. Fluid Mech., vol. 32, 2000, pp. 613–657.

8. Dommermuth, D. G.; Rottman, J. W.; Innis, G. E.; and Novikov, E. A.: Numerical

simulation of the wake of a towed sphere in a weakly stratified fluid. J. Fluid Mech.,

vol. 473, 2002, pp. 83–101.

9. Brucker, K. A.; and Sarkar, S.: A comparative study of self-propelled and towed wakes

in a stratified fluid. J. Fluid Mech., vol. 652, 2010, pp. 373–404.

10. Lin, J. T.; and Pao, Y. H.: Wakes in stratified fluids. Ann. Rev. Fluid Mech., vol. 11,

1979, pp. 317–338.

11. Godeford, F. S.; and Cambon, C.: Detailed investigation of energy transfer in homogeneous stratified turbulence. Phys. Fluids, vol. 6, no. 6, 1994, pp. 2084–2100.

12. Spedding, G. R.: Vertical structure in stratified wakes with high initial Froude number.

J. Fluid Mech., vol. 454, 2002, pp. 71–112.

13. Gourlay, M. J.; Arendt, S. C.; Fritts, D. C.; and Werne, J.: Numerical modeling of initially turbulent wakes with net momentum. Phys. Fluids, vol. 13, no. 12, 2001, pp. 3783–

3802.

14. Redford, J. A.; Castro, I. P.; and Coleman, G. N.: On the universality of turbulent

axisymmetric wakes. J. Fluid Mech., vol. 710, 2012, pp. 419–452 (referred to herein as

RCC).

42
15. Bonnier, M.; Bonneton, P.; and Eiff, O.: Far-wake of a sphere in a stably stratified fluid.

Appl. Sci. Res., vol. 59, 1998, pp. 269–281.

16. Pal, A.; de Stadler, M. B.; and Sarkar, S.: The spatial evolution of fluctuations in

a self-propelled wake compared to a patch of turbulence. Phys. Fluids, vol. 25, 2013,

p. 095106.

17. Meunier, P.; and Spedding, G. R.: A loss of memory in stratified momentum wakes.

Phys. Fluids, vol. 16, 2004, pp. 298–305.

18. de Stadler, M. B.; and Sarkar, S.: Simulation of a propelled wake with moderate excess

momentum in a stratified fluid. J. Fluid Mech., vol. 692, 2012, pp. 28–52.

19. de Stadler, M. B.; Sarkar, S.; and Brucker, K. A.: Effect of Prandtl number on a

stratified turbulent wake. Phys. Fluids, vol. 22, 2010, p. doi:10.1063/1.3478841.

20. Gottlieb, D.; and Orszag, S. A.: Numerical analysis of spectral methods: theory and

applications. SIAM, 1977.

21. Kim, J.; Moin, P.; and Moser, R.: Turbulence statistics in fully developed channel flow

at low Reynolds number. J. Fluid Mech., vol. 177, 1987, pp. 133–166.

22. Spalart, P. E.; Moser, R. D.; and Rogers, M. M.: Spectral methods for the Navier–Stokes

equations with one infinite and two periodic directions. J. Comput. Phys., vol. 96, no. 2,

1991, pp. 297–324.

23. Rogallo, R. S.: Numerical experiments in homogeneous turbulence. NASA Technical

Memorandum 81315, 1981.

24. Shariff, K.; Verzicco, R.; and Orlandi, P.: Three-dimensional vortex ring instabilities.

J. Fluid Mech., vol. 279, 1994, pp. 351–375.

25. Spedding, G. R.; Browand, F. K.; and Fincham, A. M.: Turbulence, similarity scaling

and vortex geometry in the wake of a towed sphere in a stably stratified fluid. J. Fluid

Mech., vol. 314, 1996, pp. 53–103.

26. Redford, J. A.; Llor, A.; and Watteaux, R.: Self-similarity of variable acceleration

Rayleigh-Taylor turbulent mixing layers. Proceedings of the ETMM9 Engineering Turbulence Measurement and Modelling symposium, ETMM9, Thessaloniki, Greece, June

6–8, 2012., 2012.

27. Mowbray, D. E.; and Rarity, B. S. H.: A theoretical and experimental investigation

of the phase configuration of internal waves of small amplitude in a density stratified

liquid. J. Fluid Mech., vol. 28, 1967, pp. 1–16.

28. Turner, J. S.: Buoyancy Effects in Fluids. Cambridge University Press, 1973.

29. Lighthill, M. J.: Waves in Fluids. Cambridge University Press, 1978.

30. Spedding, G. R.: Anisotropy in turbulence profiles of stratified wakes. Phys. Fluids,

vol. 13, 2001, pp. 2361–2372.

31. Townsend, A.: The Structure of Turbulent Shear Flow. Cambridge University Press, 2

ed., 1976.

32. Itsweire, E. C.; Koseff, J. R.; Brigss, D. A.; and Ferziger, J. H.: Turbulence in stratified

shear flows: implications for interpreting shear-induced mixing in the ocean. J. Phys.

Oceanogr., vol. 23, 1993, pp. 1508–1522.

43
33. Bradshaw, P.: Effects of extra rates of strain – review. Zoran Zaric Mem. Seminar,

Dubrovnik. Hemisphere, 1988.

34. Uberoi, M. S.; and Freymuth, P.: Turbulent energy balance and spectra of the axisymmetric wake. Phys. Fluids, vol. 13, no. 9, 1970, pp. 2205–2210.

35. Johansson, P. B. V.; George, W. K.; and Gourlay, M. J.: Equilibrium similarity, effects

of initial conditions and local Reynolds number on the axisymmetric wake. Phys. Fluids,

vol. 15, no. 3, 2003, pp. 603–617.

44
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

Direct Numerical Simulation of a Weakly Stratified Turbulent Wake 

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Redford, J. A.; Lund, T. S.; Coleman, Gary N.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20448

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 01

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

Direct numerical simulation (DNS) is used to investigate a time-dependent turbulent wake evolving in a stably stratified background. A large initial Froude 

number is chosen to allow the wake to become fully turbulent and axisymmetric before stratification affects the spreading rate of the mean defect. The 

uncertainty introduced by the finite sample size associated with gathering statistics from a simulation of a time-dependent flow is reduced, compared to earlier 

simulations of this flow. The DNS reveals the buoyancy-induced changes to the turbulence structure, as well as to the mean-defect history and the terms in the 

mean-momentum and turbulence-kinetic-energy budgets, that characterize the various states of this flow -namely the three-dimensional (essentially 

unstratified), non-equilibrium (or 'wake-collapse') and quasi-two-dimensional (or 'two-component') regimes observed elsewhere for wakes embedded in both 

weakly and strongly stratified backgrounds. The wake-collapse regime is not accompanied by transfer (or 'reconversion') of the potential energy of the 

turbulence to the kinetic energy of the turbulence, implying that this is not an essential feature of stratified-wake dynamics. The dependence upon Reynolds 

number of the duration of the wake-collapse period is demonstrated, and the effect of the details of the initial/near-field conditions of the wake on its 

subsequent development is examined.

15. SUBJECT TERMS

Direct numerical simulation; Stratification; Turbulence; Turbulent wakes; Wakes

18. NUMBER

 OF 

 PAGES

49

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

794072.02.07.02.03 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TM-2014-218523

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 09 - 2014