12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

THE EFFECT OF SWEEP-ANGLE VARIATION ON THE TURBULENCESTRUCTURE IN HIGH SPEED TURBULENT REACTIVE FLOWSRadouan Boukharfane
Safae El Misaoui
Mohammed VI Polytechnic University Mohammed VI Polytechnic University
MSDA group, Benguerir, Morocco
MSDA group, Benguerir, Morocco
radouan.boukharfane@um6p.ma
safae.elmisaoui@um6p.ma

Matteo Parsani
Nilanjan Chakraborty
King Abdullah University of Science and Technology School of Engineering, Newcastle University
CEMSE, ECRC, Thuwal, KSA
Newcastle-Upon-Tyne, United Kingdom
matteo.parsani@kaust.edu.sa
nilanjan.chakraborty@newcastle.ac.uk

ABSTRACT

ulation (DNS) of flow over a swept rearward-facing step and
explored the changes in the statistical turbulence structure due
to a variation of a sweep angle from to 0◦ up to 60◦ . More recently, Meldi et al. (2020) carried out a parametric variations
applied both on the sweeping angles and the shear intensity parameter applied at the inlet of an incompressible mixing layer.
Despite previous interest in supersonic combustion, few compressible reacting shear flows have been investigated recently.
On the one hand, a few experiments based on a configuration
close to a mixing layer with convective Mach numbers below
0.4 have been conducted (Cheng et al., 1994). On the other
hand, few compressible reacting shear flows have been investigated recently To address the above issues in the flow transition and self-similarity regions and address some of the inconsistent conclusions in comparison with experiments about
the Reynolds stress tensor evolution in the self-similar region,
a comprehensive DNS study of a compressible mixing layer at
Mac = 0.48 is conducted with the objective to extend the previous works of Boukharfane et al. (2021a) and Boukharfane
et al. (2021b) by assessing the relative impact of the inlet skew
angle on the turbulence structure and the flow evolution whileconsidering heat release effects, multicomponent mixtures anddetailed transport models.

The relevance of quantifying heat release and compressibility effects in turbulent compressible spatially developing
shear flows using detailed kinetics as well as detailed multicomponent transport description is of important relevance inthe operation of supersonic air-breathing engines, whose efficiency closely depends on the molecular mixing between oxidizer and fuel. The present study is devoted to investigated
how the variation of a single parameter applied at the inlet,
namely the sweeping (or skewness) angle, affects the turbulence structure and the evolution of the flow in a compressible
and multi-component inert and reacting mixing layers using
detailed chemistry and multicomponent transport. The preliminary computations indicate that the transitional region and thefully developed turbulence region depends strongly on the degree of flow sweeping at the inlet. The qualitative inspection of
the results indicates the sweeping angle tends to increase more
quickly the growth of the inlet structures, leading therefore toan enhancement the development of the mixing region. Thisis confirmed with the time-averaged statistics that characterize the mixing layer development. One possible perspective of
the present work concerns flow control of mixing layer setup
and reliable comparison between experiment, simulations andturbulence modelling.

INTRODUCTION

1

A significant amount of research has been done studying
various aspects of reactive compressible turbulent shear flows
in which the two mixing layers emerged between the injectedfuel and the surrounding supersonic airflow within the combustion chamber are assumed to be parallel. However, sweeping effects due to controllable or uncontrollable lack of perfectalignment in practical applications or even experiments, which
can affect the development of turbulence, have been addressedmuch more rarely, even though the the concept of enhancing
mixing through the introduction of a swirl effect is quite old.
Swithebank & Chigier (1969) considered several skew angles
and concluded that the sweeping effect is responsible for increased mixing in the incompressible turbulent round jet flow
and might cause similar results in compressible flows as well.
Kaltenbach (2003) performed a series of direct numerical sim-

NUMERICAL SOLVER

The unsteady, three-dimensional, compressible NavierStokes equations are considered for a multi-species reactivegas mixture:∂t ρ + ∂ j ρu j = 0 ,∂ ( ρu ) + ∂ ρu u  + ∂ p = ∂ τ ,
t
i
j
i j
i
j ij


∂
(
ρe
)
+
∂
(
ρe
+
p
)
u
=
∂
u
t
t
j
i
j
i τi j − ∂ j q j , t∂ ( ρY ) + ∂ ρu Y  = −∂ ρY u  + ρ ω̇ .
t

α

j

j

α

j

α αj

α

where

q = −λ ∂ T + ∑Nsp ρYα u ( hα + RT χeα /Wα ) ,
j
j
α=1

 αj
Xβ −YβeρYα uαi = − ∑Nsp ρ D∂
pj
αβ ∂ j Xβ +
p
β =1

1

(1a)
(1b)
(1c)
(1d)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

0

Nre

00
0
ω̇α = ∑ να,i− να,i
i=1



f

Nsp

0

Nsp

00

ki ∏ [Xα ]να,i − kri ∏ [Xα ]να,i
α=1

,

x2

α=1

This system of conservative equations is based on the following set of variables. The letter t represents the time and thesymbols ∂t and ∂i denote the partial differential operators ∂ /∂tand ∂ /∂ xi , respectively. The velocity component in direction
i is denoted by ui , the density by ρ, the pressure by p, the totalenergy (the sum of the internal specific energy, e, and the kinetic energy) by et = e + ui ui /2, and the mass fraction of the
α th species by Yα (with α = 1, · · · , Nsp , where the integer Nsp
denotes the number of species). Herein, the density, pressure,and temperature, T , are interrelated through the ideal gas law,p = ρRT /W , where W −1 = ∑N
α=1 Yα /Wα and Wα is themolecular mass of of the α th species.
The stress tensor is eval
uated as τi j = µ ∂ j ui + ∂i u j + ( κ − 2µ/3 ) ∂k uk δi j , whereµ refers to the dynamic viscosity and κ to the bulk viscosity.
λ stands for the thermal conductivity and hα for the sensible
00 and k f stand, respecenthalpy of the chemical specie α. να,iitively, the forward stoichiometric coefficient and the forward
rate constant of the i-elementary reaction (Nre being the num0 and kr are their reverse counber of reactions), whereas να,ii
terparts. The rate constants are assumed to obey an Arrhenius
law and [Xα ] are the molar concentrations of species α.
The viscous and molecular diffusion flux functions arecomputed using an eighth-order centred finite difference
scheme. The temporal integration is performed by using athird-order total variation diminishing (TVD) Runge–Kutta algorithm (Gottlieb & Shu, 1998). The enthalpies and specificheat capacities at constant pressure of each species are expressed as polynomial forms of temperature using a series ofcoefficients determined from the JANAF databases (Prophet
et al., 1971). A detailed verification of the solver may befound in Boukharfane (2018), which assembles five elementary and verification subsets including the Sod’s shock tube
problem, the ignition sequence of a multi-species mixture in a
shock tube, a one-dimensional laminar premixed flame, isentropic vortex, and shock–vortex interactions.

2

20
=5

!

U

eed

x1
ter

it
Spl

p

L x1

x2

h

h-sp
Hig

ream

airst

x1
ing

Mix

ζ

late

× δω

x3

r
laye

am

e
irstr

ed a

spe
ow-

L x2 = 120 × δω0

and

L
l

U

Figure 1: Schematic of the computational arrangement

often identified through the convective Mach Number, Mac =∆U/(ch + cl ). For the present study, it is set to Mac = 0.48.
A schematic of the computational arrangement is shown in
Fig. 1. The size of the computational domain for each caseis set as (Lx1 × Lx2 × Lx3 ) /δω0 = 520 × 120 × 60 resolved by
Nx1 × Nx2 × Nx3 = 2300 × 400 × 200 grid points. In addition
to the reference case, corresponding to ζ = 0◦ , three numerical simulations of mixing layers are conducted corresponding
to three angles of sweeping, ζ = 5◦ , 10◦ and 15◦ .

2.1
2.2

NUMERICAL RESULTS
MIXING LAYER FLOW FIELD CHARACTERISTICS

The classical shapes and organization of the large-scale
turbulent structures in the mixing layer is reported in fig. 6aand section 2.5. Compressibility effects are visible in thepresent setup due to absence of spanwise Brown-Roshko
rollers and the presence of the three-dimensional and disorganized turbulent structures with a dominance of small-scale eddies. Heat release has also a noticeable impact in the transitionof the mixing layers. Indeed, complex oblique flow structures
grow and the formation of two-dimensional large structures isreadily observed through the streamwise direction as a consequence of the relaminarizing effect induced by heat release;see fig. 6a and section 2.5. On the other hand, the sweeping mixing layers featuring a sweeping angle of ζ = 15◦ show
a different flow pattern compared to the case ζ = 0. In thefirst part of the mixing process, which features mainly twodimensional spanwise vortices, the axis of the coherent structures in the case ζ = 15◦ is titled with respect to the normal
incoming flow direction. In the unskewed case, the shear layers roll up into Λ-shaped vortices, which quickly evolve into
hairpin vortices and then sheaths with intense vorticity andthen break up into small slender vortices. The skewed case
exhibits visibly no Λ-shaped vortices, only elongated sheaths
are present near the inlet of the mixing layers that lately breakup. In the second region characterized by the break-downof two-dimensional primary coherent structures to small-scale
structures, the vortex roll-up process takes place early and thegrowth of inlet structures are more magnified in the presence ofsweeping effects. Substantial increase of mixing-layer growthwith the presence of skewness effects is observed. Moreover,the interactions between large-scale rollers and rib vortices
that trigger the turbulence transition and the break down tosmall-scale structures occurs earlier with ζ = 15◦ comparedto ζ = 0◦ .

SIMULATION SETUP

The turbulent spatially evolving mixing layers is a flow
originating from the interaction between two homogeneous
streams having different asymptotic, high Uh (or hot oxidizer
stream O2 /H2 O/N2 ) and low Ul (or cold fuel stream H2 /N2 )speed stream velocities. The sweeping of the spatially mixing layers indicate a lack of perfect alignment between the twoparallel and asymptotic velocities by a constant angle ζ mea-
sured on the x1 –x3 planes in such manner that u h = Uh , 0, 0and u l = Ul cos ζ , 0, Ul sin ζ . The flow at the inlet (x1 = 0)
is a hyperbolic tangent profile for the streamwise and spanwise velocities, species mass fractions and density/ All numerical simulations are initialized in the computational domainusing the unperturbed hyperbolic tangent profile for the velocity, species mass fraction and density. The Reynolds number based on the initial vorticity thickness, the velocity difference ∆U = Uh − Ul , and the average dynamic viscosity
is 640. Dirichlet boundary conditions are applied at the two
supersonic inlets. Perfectly non–reflecting boundary conditions are set at the outflow and the cross–stream directions,whereas periodic boundary conditions are settled along the
spanwise directions. For equal densities and specific heats,
the effects of compressibility on the shear-layer growth rate is2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022(a) ζ = 00◦(b) ζ = 15◦
Figure 2: Q-criterion colored with local Mach number to visualize the development of coherent structures in the mixing

layer for both angles of skewness.
2.3

GROWTH RATES

To have a rough estimation of the local size of the turbulent structures in the mixing layer, in Fig. 3 we show the momentum and vorticity thicknesses, normalized by ηδω0 , with
η = u h −uul / u h +uul . Here, the momentum thickness isdefined as follows!
ZZ
hu1 i f − Uhhu1 i f − Uh
1 +∞ +∞
hρi h
1−dx2 dx3 ,δθ =ρ0 −∞ −∞
u −uulu h −uul(2)
where ρ0 stands for the arithmetic mean value between thedensities of the oxidizer and fuel inlet streams ρ0 = (ρ h +ρ l )/2. The mixing layer evolves to a self-similar state past
x1 /δω0 ≈ 200, regardless of the skew angle, whence near linear growth is observed. The vorticity thickness shows significant differences in the transition region, in which skewing promotes faster growth. Heat release has also a noticeable impact
in the transition of the mixing layer since the reacting cases
shows higher growth compared to the inert case as a consequence of the relaminarizing effect induced by heat release.
The slopes of the vorticity thickness regression lines indicate
small but consistent growth with the skew angle is observed
also in the fully developed region, which is likely to be exploitable in Scramjet engines applications. The momentum
thickness exhibits much less dependence on ζ , with exceptionof the the highest value under consideration, ζ = 15◦ and heatrelease tend to noticeably reduce the growth rate in the fully
developed region of mixing layers.

2.4(a) ζ = 00◦

REYNOLDS STRESSES

Further insight into the effects of flow skewing may begained from inspection of distribution of the Reynolds stresses,
which are also expected to reach an asymptotic self-similar
regime. The streamwise and spanwise Reynolds stress profiles are shown as a function of the self-similar cross-stream
coordinate x2 /δω in Fig. 4. Both stresses exhibit tendency
to collapse, and the streamwise stress profiles show moderate asymmetry, peaking in the region occupied by the heavier
oxidizer (lower density region), as noted by Brown & Roshko(1974) and Pantano & Sarkar (2002). Although the unskewed

Figure 3: Comparison the streamwise evolution of (a)vorticity thickness and (b) momentum thickness of themixing layer. The dashed lines denote fitted linear regression curve.

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022(a) ζ = 00◦(b) ζ = 15◦(c) ζ = 00◦(d) ζ = 15◦

wise location for all the values of the components of Reynold
stress, as in the present compressible mixing layers. The
skewed case exhibits a peak values for streamwise turbulence
intensity that start to rise upstream further compared to theunskewed case, but a lower peak value is reached. However,the streamwise intensity relaxes monotonically to an asymptotic value similar to the unskewed case, meaning that the peakvalue of R11 seems to be unaffected by skewing. However, thespanwise turbulent intensity experiences a sharp increase anddecays to an asymptotic higher than R11 . A possible explanation for these findings is that as the skew angle is increased, theelongation of large-scale turbulence structures become higher
in the x1 -direction relative to x2 and x3 directions, which tends
to be compensated further downstream by the increase of thefluctuations of R22 and, a mainly of R33 . Beyond the maximum peaks, the turbulent intensities of R22 and R33 seemto be suppressed in the absence of skewing effect much faster,which indicates that skewed mixing layers induce more intense
turbulence and enlarge the the momentum exchange between
the turbulent structures and the inflow mainstream.In the mean
convection direction, the components involving of spanwise
direction x3 present lower turbulent intensity peaks comparedto the natural direction for the most skewed case.

Figure 4: Profiles of streamwise (top row) and spanwise(bottom row) Reynolds stresses at various streamwise
locations for ζ = 0◦ and 15◦ . Solid lines corresponds tonatural direction, while dashed lines stands for the mean
convective direction.

case exhibits higher levels of the streamwise Reynolds stress inthe development regime, skewed and unskewed cases achieveself-similarity at approximately the same streamwise location,
x1 /δω0 = 450, and the maximum√ level in the self-similar stateis similar. The peak value of R11 / uh − ul in the present
simulation of 0.137 is close to the numerical data of Debisschop et al. (1994), who reported a value 0.141 at convective Mach number Mac = 0.525. Stronger √effect of skewing
is found on the spanwise Reynolds stress R11 / uh − ul ,for which the skewing yields substantially higher peak values both in the development and√in the self-similar region. In
the mean convection direction, R33 / uh − ul exhibits nat√urally lower amplitude, whereas R11 / uh − ul seems to benearly independent of the system of coordinate in where it iscomputed.
The spatial evolution of the peak values of turbulent intensity distributions along the normalized streamwise x1 -direction
is displayed in Fig. 5. The profile of peak intensities for both

2.5

ANISOTROPY EFFECT

A key tool in the development of advanced turbulence
closures is the consideration of the anisotropy effects (Vreman et al., 1996), which can be reformulated using either the
Reynolds stress anisotropybK
ij =Ri j − 23 δi j K
,
2K

(3)

or the dissipation anisotropy as
bεi j =

εi j − 23 δi j ε
,
2ε

(4)

where δi j is the Kronecker delta. The streamwise evolutionof anisotropy tensor, where Ri j and K are integrated along
the spatial mixing layer are displayed in Fig. 6 for ζ = 00◦
and ζ = 15◦ . The similarity between bK and bε (not shownhere) profiles indicates that the flow organization fairly influences the small scale turbulence structure, which induces almost the same anisotropy in both quantities. The other cases
are not shown since the trends observed tends to amplify as ζ
becomes bigger. The profiles of anisotropies is quite similar
either using the Reynolds stress or dissipation tensors. In the
pre-transitional regime, all the anisotropies increase to reacha peak value before converging to an almost asymptotic constant value in the self-similar region. Skewing the inlet inflow
strongly affects the initial evolution of the anisotropies and alK
ters their asymptotic values. As ζ increases, the bK
33 and b23
profile increase above zero, indicating that a higher level ofskewing induces an increase of the level of anisotropy. The
use of the isotropic assumptions will likely be inaccurate at
higher skew angle ζ .

Figure 5: Peak values of turbulent intensity distributions

along the normalized x1 -axis in (a) the natural directionand (b) in the mean convective direction. Solid lines:
ζ = 00◦ , and dashed lines: ζ = 15◦ .

Given that the tensor bK has zero trace, Lumley & Newman (1977) showed that all the anisotropy states are describedby the two invariants of bK given byskewed and unskewed cases has a distribution different fromthe classical incompressible mixing layers, in which the peakintensities are achieved at almost the same normalized stream-

IIK = −4

1 2
Λ1 + Λ22 + Λ23 ,2(5)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022(a) ζ = 00◦(b) ζ = 05◦(c) ζ = 10◦(d) ζ = 15◦(a) ζ = 00◦

Figure 7: Reynolds stress anisotropy invariants using

Lumley triangle constraints colored by the streamwise
coordinate.
angle (axisymmetric disk-like state or oblate spheroid shape)as the skew angle becomes higher. These findings can also
be observed by considering the determinant of the normalized
Reynolds stress tensor F, which is related to IIK and IIIK as
follows
F = 1 + 9IIK + 27IIIK .

with F being one denote a perfectly isotropic and threedimensional turbulence, and F being zero denote a pure twocomponent turbulence and completely anisotropic. As shownin Fig. 8, all cases almost collapse on the self-similarity region (giving another statistical tool to measure self-similarity)
and are close to unity, indicating an approach towards isotropy.
None of the cases reach unity, meaning that the present flows
exhibit anisotropy throughout the fields. The main effect ofskewing occurs between the inlet of the computational domain and the transition to turbulence. In particular, the level isisotropy is lower as the skew angle becomes higher.

Figure 6: Longitudinal evolution of the anisotropy ten-sors b K .

and
IIIK =

1 3
Λ1 + Λ32 + Λ33 ,
3

(7)

(6)

where Λ1 , Λ2 , and Λ3 are the eigenvalues of bK . Physically, IIK is an indicator of anisotropy in the flow, while
IIIK denotes the nature of anisotropy. The different realizableanisotropic states of the turbulent velocity flowfield are contained within a triangle in the IIK -IIIK plan, the so-called
Lumley triangle Lumley (1979). One should caution that the
topology of of the energy distribution induced from the analysis of the IIK -IIIK plan should not be interpreted as theshape of the flow eddies since the characteristic shapes associated to the stresses and and those to the eddies are not generally
equivalent (Simonsen & Krogstad, 2005). The resulting invariant mappings along the streamwise direction are portrayed inFig. 7 for the skewing angles with the Lumley triangle boundaries labelled. The limiting states are shown by the threesolid lines. In all the cases, the turbulence is essentially twocomponent in the transitional regime and three-component at
the self-similarity state for all the case. The main difference
occurs at the pre-transitional regime. Indeed, for the unskewed
case and as the streamwise coordinate increases, the turbulenceanisotropy is relatively approaching the axisymmetric expansion bound (rod-like shapes) and the isotropic state. However, as ζ becomes higher, the invariant mapping get closer
to the one-component limit along the axisymmetric bound to
the right. The flow slightly tends to be in close proximity tothe plane-strain condition to the left-side of the Lumley tri-

Figure 8: Longitudinal evolution of the invariant function

F.

3

CONCLUSION

In the present paper, structural evolution and turbulent
statistics of a compressible mixing layer with Mac = 0.48 areinvestigated employing direct numerical simulation to identify

512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
of homogeneous turbulence. Journal of Fluid Mechanics
82 (1), 161–178.
Meldi, M., Mariotti, A., Salvetti, M. V. & Sagaut, P. 2020 Numerical investigation of skewed spatially evolving mixing
layers. Journal of Fluid Mechanics 897.
Pantano, C. & Sarkar, S. 2002 A study of compressibility effects in the high-speed turbulent shear layer using direct
simulation. Journal of Fluid Mechanics 451, 329.
Prophet, H., Stull, D. R., United States & National Bureau ofStandards 1971 JANAF thermochemical tables. U.S. Dept.
of Commerce, National Bureau of Standards.
Simonsen, A. J. & Krogstad, P.-Å. 2005 Turbulent stressinvariant analysis: Clarification of existing terminology.
Physics of Fluids 17 (8), 088103.
Swithebank, J. & Chigier, N. A. 1969 Vortex mixing for supersonic combustion. In Symposium (International) on Combustion, , vol. 12, pp. 1153–1162. Elsevier.
Vreman, A. W., Sandham, N. D. & Luo, K. H. 1996 Compressible mixing layer growth rate and turbulence characteristics.
Journal of Fluid Mechanics 320, 235–258.

the effect of skewness. Both instantaneous and time-averaged
data are utilized to better understand the flow dynamical behaviors in the presence or not of skewness. The instantaneous
visualization results indicate an amplification of the inlet disturbances resulting in a faster increase of the mixing process.
The higher mixing layer growth rate is a direct consequenceincreaed velocity fluctuations, i.e., increased turbulent kinetic
energy and increased Reynolds stresses. Their increase withskewing effect is traced back to a decrease of the pressurestrain correlations and an increase of the production rate of the
spanwise Reynolds stress. These correlations mainly steer the
redistribution mechanism of turbulent fluctuating energy fromthe streamwise Reynolds stress component to the spanwise
component. The peak streamwise normal stress showed little
variation, while the transverse normal, spanwise normal, andprimary shear stresses all increased in peak magnitude with increasing skew angle. In addition, Reynolds stress anisotropy
observations were made. From the similarity profiles, the main
effect of skewing were observed in the transition region. The
invariant mappings of the mixing layers were found to tend toward the one component limit of the Lumley triangle as the
skew angle is increased. A slightly reduction of all main terms
in the Reynolds stress transport equations and the transport
equation of the turbulent kinetic energy with increasing the
skew angle ζ has been noticed, and an analysis of the corresponding transversely integrated terms has shown that thepressure-strain correlation of the streamwise Reynolds stressis directly connected to the momentum thickness growth rate
and responsible for its increase.

REFERENCES
Boukharfane, R. 2018 Contribution à la simulation numériqued’écoulements turbulents compressibles canoniques. PhD
thesis, École Nationale Supérieure de Mécanique etd’Aérotechnique.
Boukharfane, R., Er-raiy, A., Elkarii, M. & Parsani, M.
2021a A direct numerical simulation study of skewed threedimensional spatially evolving compressible mixing layer.
Physics of Fluids 33 (11), 115114.
Boukharfane, R., Er-Raiy, A., Parsani, M. & Hadri, B. 2021b
Skewness effects on the turbulence structure in a high-speed
compressible and multi-component inert mixing layers. In
AIAA Aviation 2021 Forum, p. 2915.
Brown, G. L. & Roshko, A. 1974 On density effects and large
structure in turbulent mixing layers. Journal of Fluid Mechanics 64 (4), 775–816.
Cheng, T. S., Wehrmeyer, J. A., Pitz, R. W., Jarrett Jr, O. &Northam, G. B. 1994 Raman measurement of mixing andfinite-rate chemistry in a supersonic Hydrogen-air diffusion
flame. Combustion and Flame 99 (1), 157–173.
Debisschop, J. R., Chambers, O. & Bonnet, J. P. 1994 Velocity field characteristics in supersonic mixing layers. Experimental thermal and fluid science 9 (2), 147–155.
Gottlieb, S. & Shu, C.-W. 1998 Total variation diminishing Runge–Kutta schemes. Mathematics of Computation
67 (221), 73–85.
Kaltenbach, H.-J. 2003 The effect of sweep-angle variation on the turbulence structure in a separated, threedimensional flow. Theoretical and Computational Fluid
Dynamics 16 (3), 187–210.
Lumley, J. L. 1979 Computational modeling of turbulent
flows. In Advances in Applied Mechanics, , vol. 18, pp.
123–176. Elsevier.
Lumley, J. L. & Newman, G. R. 1977 The return to isotropy

6