Elementos de Turbulencia Etc.

Download as pdf or txt
Download as pdf or txt
You are on page 1of 88

Università degli Studi della Campania Luigi Vanvitelli

Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

ELEMENTS OF TURBULENCE

ALDO TAMBURRINO
Associate Professor Universidad de Chile
Visiting Professor Università della Campania Luigi Vanvitelli

INTRODUCTION
The motion of a fluid can be characterised by different regimes, depending
on the specific feature that we are interested in. Thus, if we are interested in the
flow dependency on the time, it could be steady (the flow does not change with
time) or unsteady (the flow depends on the time). If we are interested in the
changes in the space, it could be uniform (velocity and flow section do not
changes in space), or spatially varied. In the case of spatially varied flow, when
the variation is small, the flow is gradually varied. In this case, the curvature of
the streamlines is small. On the contrary, when the streamlines have strong
curvature, the flow is named rapidly varied. A property of the gradually varied
flows is that the pressure distribution can be considered as hydrostatic. There
are many other features of the flow that we could want to study, and we could
define another regimes. However, we will focus in the motion of the fluid
particles. That motion can be ordered, the flow moving in “layers”, and the flow
regime is called laminar. On the contrary, the motion of the particles can be
disorganised or random and the regime is named turbulent. Of course, there is
not an abrupt change from the laminar to the turbulent regime, but it goes
through a transition laminar-turbulent.
Osborne Reynolds published in 1883 an article in which he described the
motion of water in a tube in which he injected dye. The experimental facility
used by Reynolds and depicted in the article is shown in Fig. 1. As a curiosity,
Fig. 2 shows the original apparatus used by Reynolds, currently in the
University of Manchester.
For low water velocities, the dye injected in the upstream end of the tube
followed a well-defined linear trajectory, as shown in Fig.3 (taken from Reynold’s
paper). Reynolds called this motion “direct”. For higher velocities, after some
distance from the entrance of the tube, the dye “at once mix up with the
surrounding water, and fill the rest of the tube with a mass of coloured water”,as
shown in Fig. 4. But the most interesting feature was observed when the tube
:: 1 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

was illuminated with an “electric spark”: the dye showed “more or less distinct
curls, showing eddies”, as those sketched in Fig.5. Because of this eddy pattern,
Reynolds called this motion also “sinuous”. In the current language we call

Fig. 1.- Experimental set-up of Reynold´s experiment, taken from his


paper from 1883.
“laminar” and “turbulent” regimes what Reynolds called “direct” and “sinuous”
motions. He also found the limits at which the regime ceased to be laminar and
become turbulent was proportional to the tube diameter 𝐷, the mean flow
velocity 𝑈, and inversely proportional to the kinematic viscosity  , that is to
say 𝑈𝐷⁄. This dimensionless parameter is known as the Reynolds number:
𝑈𝐷
𝑅𝑒 = (1)

The accepted values of the Reynolds numbers for the limits of the different
regimes are:

:: 2 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

- Laminar regime : 𝑅𝑒 < 2000


- Laminar-Turbulent transition regime : 2000 < 𝑅𝑒 < 4000
- Turbulent regime : 𝑅𝑒 > 4000

Presence of eddies in the flow was not something new when Reynolds
performed his experiment. The eddy motions in water flows were already

Fig. 2.- The original experimental facility used by Reynolds, now


at the University of Manchester (Picture taken in May, 2015)

depicted in a series of drawings by Leonardo da Vinci, as it is shown in Fig. 6,


who in the description of the motion of the surface of the water stated that it has

:: 3 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

two components, one given by the “impetus” of the main motion and another
associated to the eddies. Leonardo was the first one that used the term
turbolenza in about 1550.
Existence of eddies in the flow has an important effect in the flow
resistance. Before than Reynolds’ experiment was already known that for low
velocities (i.e., laminar regime), the flow resistance was proportional to the mean
velocity (or discharge), and for higher velocities (i.e., turbulent regime) it was
proportional to the square of the mean velocity.

Fig. 3.- Flow pattern of the dye at low velocities. “Direct”


motion or laminar regime.

Fig. 4.- Flow pattern of the dye at high velocities.


“Sinuous” motion or turbulent regime.

Fig. 5.- Flow pattern of the dye in the “sinuous” motion or


turbulent regime when illuminated with a sparkling light.
Reynolds observed the eddy or curly motion of the fluid.

WHAT IS TURBULENCE?
It is difficult to answer this question. As the velocity is continuously
changing (o fluctuating) in time, we can say that it is an unsteady phenomenon.
However, usually we are not interested in the instantaneous values of the
velocity (or pressure) of a turbulent flow, but we want to know some time
averaged values. Thus, the continuous variation in time of the flow properties

:: 4 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

has originated a statistical description of the flow in terms of temporal mean


values.

Fig. 6.- Drawing of Leonardo da Vinci showing the presence of


eddies in a flow with free Surface.

Due to the complexity of the turbulent flows, most of what we know about
them comes from experiments. By means of experiments, using visualization
techniques or direct measurements of the variables that we are interested in
(usually, velocity and pressure), we can get a general description of the processes
involved in the turbulent flows. The use of the experimental approach does not
mean that theories have not been developed, or that are of lesser value than the
experiments. On the contrary, the theoretical analysis guide us regarding the
variables that need to be measured for a better understanding of the turbulence
phenomena. Obviously, as in all branches of science, the interaction between
theory and experiments is necessary, and both of them help to increase our
understanding of the turbulence. It is important to remark that turbulence still
is an open problem. That is to say that currently the problem has not been solved
starting from the principles of the physics. Thus, any theory on turbulence will
always rely on some experimental data.
Up to now, there is not a definition of turbulence that can describe it
completely. Von Karman (1937) defined it as “an irregular motion which in
general makes its appearance in fluids, gaseous or liquid, when they flow past
solid surfaces or even when neighbouring streams of the same fluid flow past or
over one another”. According to Hinze (1959) “turbulent fluid motion is an
irregular condition of the flow in which the various quantities show a random
variation with time and space coordinates, so that statistically distinct average
values can be discerned”. This definition has a couple of consequences: the

:: 5 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

variation must be space-temporal, and it can be studied by means statistical


analysis.
We can look for the definition given by different dictionaries, for example:
- English
o Oxford Dictionary (https://en.oxforddictionaries.com/): “Violent or
unsteady movement of air or water, or of some other fluid”
o Dictionary.com (http://www.dictionary.com): “The haphazard
secondary motion caused by eddies within a moving fluid”
o Merriam-Webster Dictionary: “Departure in a fluid from a smooth
flow”

- Italian
o Dizionario Italiano (www.dizionario-italiano.it): “Fenomeno per cui
in un fluido liquido o gassoso in moto, in determinate condizioni, il
moto delle particelle elementari cessa di essere regolare e diventa
soggetto a forti fluttuazioni della velocità e a moti vorticosi e a
mancanza di regolarità nella traiettoria delle particelle”
o Grande Dizionario Italiano di Gabrielli Aldo
(www.grandidizionari.it): “Moto irregolare generalmente rilevabile
nei fluidi”
o Dizionari di Italiano Corriere della Sera
(http://dizionari.corriere.it/dizionario_italiano): “Moto disordinato
di un fluido, con formazione di vortici”

- Spanish
o Diccionario de la Real Academia Española (www.rae.es): “Zona en
que se desarrolla un movimiento turbulento”.
o Turbulento: “Dicho del movimiento: Propio de un fluido en el
que la presión y la velocidad fluctúan muy irregularmente en
cada punto, con la consiguiente formación de remolinos.”
As it can be seen, even the definition of the word turbulence in different
languages and dictionaries is not unique. However, from in the definitions
presented before we can group words that transmit more or less the same
meaning. For example, we find:
- “Violent or unsteady movement”, “forti fluttuazioni della velocità”,
- “Il moto delle particelle elementari cessa di essere regolare”, “Moto
irregolare”, “Moto disordinato”, “la presión y la velocidad fluctúan muy
irregularmente”

:: 6 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

- “Eddies”, “moti vorticosi”, “Moto …con formazione di vortici”,


“movimiento …con … formación de remolinos”
We note that the concept of eddies, vortici, remolinos, appears in several
definitions. As we will see later, this is a very important concept in turbulence,
used to model turbulent flows.
Given the difficulties to define turbulence, it is easier to give some
characteristics of turbulent motion. Thus, according Tennekes and Lumley
(1972), the following characteristics can be mentioned:
Turbulence characterises a kind of flow: Turbulence is not a property of a fluid,
it is a kind of flow. For this reason, the most important characteristics of
turbulent flows are not controlled by the molecular properties of the flow.
Large Reynolds numbers: Turbulent flows occur at high Reynolds numbers.
Turbulence often is originated by the instability of a laminar flow when the
Reynolds number surpasses a critical value. In turbulent flows, the effect of the
inertia of the fluid motion overcomes the stabilizing effect of the viscous forces.
Continuum: Turbulence is a continuum phenomenon, governed by the equations
of fluid mechanics. The smallest length scales that can arise in a turbulent flow
are much larger than any molecular length scale.
Irregularity: As it is observed in Reynolds’ experiment, the trajectory of a fluid
particle is highly irregular, making practically impossible to describe it by means
of deterministic analysis, and statistical methods are used.
Diffusivity: Turbulent flows are highly diffusive. For this reason, exchange rates
of mass, momentum and energy are much higher than those associated to
laminar flows. From the applications point of view, this is the most important
characteristic of turbulent flows.
Dissipation: Turbulent flows always dissipate energy. The work performed by
the viscous stresses when deform a fluid element increases the internal energy
of the fluid in detriment of the turbulent kinetic energy. Turbulence requires a
continuous supply of energy in order to sustain the viscous losses. If energy is
not supplied, the turbulence decays.
Tri-dimensional fluctuations of vorticity: Turbulence is rotational and tri-
dimensional. It is characterized by large vorticity fluctuations. An important
mechanism for the maintenance of the vorticity is “vortex stretching”. Turbulent
fluctuations of vorticity cannot be sustained by themselves in two-dimensional
flows.

:: 7 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

REVIEW OF THE BASIC EQUATIONS OF FLUID MECHANICS


General speaking, there two approaches to analyse fluid flows: the
integral approach and the differential approach. In the integral approach we use
a control volume to characterise the fluid flow (or the effect of the fluid flow) in
a finite region of the space. Some gross quantities are obtained when this
approach is used: mean velocity trough a section, force acting on a surface, total
energy exchange, etc. On the contrary, in the differential approach the
information that we get is in a point within the flow domain.
Up to now, the main principles of the physics applied to the analysis of the
fluid flows using the integral approach. Those principles are: conservation of
matter, Newton´s second law, and the first principle of thermodynamics.
Without entering in details, the resulting equations that describes those
principles, when applied to an incompressible fluid, are:
Continuity equation: It describes mathematically the law of conservation of
the matter. As the amount of matter remains constant (this principle is valid in
classical mechanics), for a fluid wit constant density, the volume rate of fluid
that enters to the control volume minus the volume rate that exits is equal to
the variation of volume of fluid in the control volume, i.e.:

𝜕𝑉
= 𝑄𝑖 − 𝑄𝑜 (1)
𝜕𝑡
In Eq. 1, 𝑉 represents the volume of the fluid inside the control volume, 𝑡
is the time, 𝑄 is the volume rate or discharge, the sub-index 𝑖 indicates the input
to the control volume and the sub-index 𝑜 the output.
If the discharge 𝑄 goes through a section with an area normal to the flow
equal to 𝐴, the average velocity is defined as

𝑄
𝑣= (2)
𝐴

In a steady state flow there is not variation in time, 𝜕𝑉⁄𝜕𝑡 = 0, and Eq. 1
is simplified to:

𝑄𝑖 = 𝑄𝑜 (3)

Eq. 3 indicates that the flow that enters to the control volume is equal to
the output flow.

:: 8 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑉𝐶

Fig. 7.- Flow through an orifice

The classical example where Eq. 1 is applied corresponds to the flow through
and orifice in a tank. Referred to Fig. 7, the question is to determine the variation
of the water depth (ℎ) in the tank in function of the time, if initially the flow
depth is ℎ0 . The transverse area of the tank is 𝐴𝑇 and the area of the orifice is
𝐴𝑜 . This simple example is repeated here to stress two isues:
i) The volume of fluid not necessarily is equal to the control volume,𝑉𝐶 ,
depicted as the segmented line in Fig. 7.
ii) Eq. 1 by itself is not enough to solve the problem.
The volume of fluid (𝑉) in the control volume 𝑉𝐶 is 𝑉 = ℎ𝐴𝑇 . There is not
flow rate entering to 𝑉𝐶 , so 𝑄𝑖 = 0, and the flow rate exiting the control volume
is 𝑄𝑜 = 𝑣𝐴𝑜 . Thus, Eq. 1 can be written as:

𝜕
(ℎ𝐴𝑇 ) = −𝑣𝐴𝑜 (4)
𝜕𝑡

We cannot go further with Eq. 4 without additional information. Up to


now, we are having two unknowns: ℎ and 𝑣, but only one equation (Eq. 4). We
need another equation involving ℎ and 𝑣. That equation will be provided by the
energy equation. But it is possible to get a result if we remember the
experimental work by Torricelli, who obtained:

𝑣 = √2𝑔ℎ (5)

Replacing Eq. 5 in Eq. 4, we obtain:

:: 9 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕ℎ 𝐴𝑜
= − √2𝑔ℎ (6)
𝜕𝑡 𝐴𝑇

Eq. 6 can be easily integrated. The integration constant is obtained from


the initial condition 𝑡 = 0, ℎ = ℎ0 .

Energy general equation: It is obtained from the first principle of the


thermodynamics, which states that the variation of energy in a system of
particles, ∆𝐸, is equal to the heat supplied to the system, ∆𝑄̂ , minus the
mechanical work done by the system, ∆𝑊. The rate of variation of energy can be
written as:

∆𝑊
∆𝐸

∆𝑄̂

Fig. 8.- First principle of thermodynamics: Variation of energy

𝑑𝐸 𝑑𝑄̂ 𝑑𝑊
= − (7)
𝑑𝑡 𝑑𝑡 𝑑𝑡

The energy of the fluids particles is given by three components:

𝐸 = 𝑈 + 𝐸𝑃 + 𝐸𝐶 (8)

where 𝑈 is the internal energy, 𝐸𝑃 is the potential energy and 𝐸𝐶 is the kinetic
energy. They are expressed as:
1
𝑈 = 𝑚𝑢 𝐸𝑃 = 𝑚𝑔𝑧 𝐸𝐶 = 𝑚𝑣 2 (9)
2
where 𝑚 is the mass, 𝑢 is the specific internal energy, 𝑔 is the acceleration due
to gravity, 𝑧 is the vertical distance from an arbitrary reference level, and 𝑣 the
:: 10 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

velocity. The specific energy (defined as energy per unit mass, 𝑒 = 𝐸 ⁄𝑚), of the
system is:
𝑣2
𝑒 = 𝑢 + 𝑔𝑧 + (10)
2
It can be shown that the variation of energy of the system can be expressed
as the change of energy per unit time of the fluid in the control volume, plus the
flow of energy through the surfaces𝑆𝐶 that define 𝑉𝐶 . Mathematically, it is
written as:
𝜕 𝑑𝑄̂ 𝑑𝑊
∫ 𝜌𝑒𝑑𝑉 + ∫ 𝜌𝑒𝑣⃗ ∙ 𝑛̂𝑑𝑆 = − (11)
𝜕𝑡 𝑉𝐶 𝑆𝐶 𝑑𝑡 𝑑𝑡
The work 𝑊 can be divided in two main components: 1) the work that the
fluid has to do when it flows (it has to overcome the forces arising from the
pressure and shear stresses), and 2) the external work (also called “shaft work”).
The external work is that done to move the shaft of a turbine (work towards the
exterior of the system, work done by the fluid, positive), or that supplied by a
pump (work towards the interior of the system, work done on the fluid,
negative). Thus:

𝑊 = 𝑊𝐸 + 𝑊𝑝 + 𝑊𝜏 (12)

where the sub-indices 𝐸, 𝑝 and 𝜏 stand for “external”, “pressure” and “shear
stresses”, respectively.
It is possible to work with the term associated to the pressure:

𝑊𝑝 = ∫ 𝑑𝑊𝑝 = ∫ 𝐹⃗𝑝 ∙ 𝑑𝑟⃗ = ∫ 𝑝𝑑𝑆𝑛̂ ∙ 𝑑𝑟⃗ (13)

𝑑𝑊𝑝 𝑑𝑟⃗
= ∫ 𝑝𝑑𝑆𝑛̂ ∙ = ∫ 𝑝𝑣⃗ ∙ 𝑛̂𝑑𝑆 (14)
𝑑𝑡 𝑑𝑡 𝑆𝐶

The work associated to the shear stresses demands to know an expression


for the shear stress, and will work on it later. Thus, Eq. 11 becomes:
𝜕 𝑑𝑄̂ 𝑑𝑊𝐸 𝑑𝑊𝜏
∫ 𝜌𝑒𝑑𝑉 + ∫ 𝜌𝑒𝑣⃗ ∙ 𝑛̂𝑑𝑆 = − − ∫ 𝑝𝑣⃗ ∙ 𝑛̂𝑑𝑆 − (15)
𝜕𝑡 𝑉𝐶 𝑆𝐶 𝑑𝑡 𝑑𝑡 𝑆𝐶 𝑑𝑡

Introducing the expression for the specific energy 𝑒 in the second term of the left
hand side of Eq. 15:

:: 11 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕 𝑣2 𝑑𝑄̂ 𝑑𝑊𝐸 𝑑𝑊𝜏


∫ 𝜌𝑒𝑑𝑉 + ∫ 𝜌 (𝑢 + 𝑔𝑧 + ) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 = − − ∫ 𝑝𝑣⃗ ∙ 𝑛̂𝑑𝑆 − (16)
𝜕𝑡 𝑉𝐶 𝑆𝐶 2 𝑑𝑡 𝑑𝑡 𝑆𝐶 𝑑𝑡

Joining the terms that are integrated over 𝑆𝐶 :


𝜕 𝑝 𝑣2 𝑑𝑄̂ 𝑑𝑊𝐸 𝑑𝑊𝜏
∫ 𝜌𝑒𝑑𝑉 + ∫ 𝜌 (𝑢 + 𝑔𝑧 + + ) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 = − − (17)
𝜕𝑡 𝑉𝐶 𝑆𝐶 𝜌 2 𝑑𝑡 𝑑𝑡 𝑑𝑡

We can see that the sum of Bernuilli was formed in the second term of Eq.
17, and it can be written as:
𝜕 𝑢 𝑑𝑄̂ 𝑑𝑊𝐸 𝑑𝑊𝜏
∫ 𝜌𝑒𝑑𝑉 + ∫ 𝜌𝑔 ( + 𝐵) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 = − − (18)
𝜕𝑡 𝑉𝐶 𝑆𝐶 𝑔 𝑑𝑡 𝑑𝑡 𝑑𝑡

where
𝑝 𝑣2
𝐵=𝑧+ + (19)
𝜌𝑔 2𝑔

Eq. 17 (or Eq. 18) is called the “energy general equation”. It can be applied
to any flow regime. Its only limitation is given by the potential energy should
derive from a gravitational field (Eq. 9).
Let’s simplify Eq. 18 assuming steady flow and considering a stream tube
as control volume, as shown in Fig. 9. The control volume has three surfaces:
section 1 (entrance), section 2 (exit) and the mantle, with surface areas 𝑆1, 𝑆2 and
𝑆𝑚 , respectively. Thus, the surface 𝑆𝐶 that defines the control volume 𝑉𝐶 can be
written as 𝑆𝐶 = 𝑆1 + 𝑆2 + 𝑆𝑚 . As the mantle is tangent to the velocity vectors,
there is not flow through its surface, and the integral over the surface of Eq. 18
is reduced to:
𝑢 𝑢 𝑢
∫ 𝜌𝑔 ( + 𝐵) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 = ∫ 𝜌𝑔 ( + 𝐵) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 + ∫ 𝜌𝑔 ( + 𝐵) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 (20)
𝑆𝐶 𝑔 𝑆1 𝑔 𝑆2 𝑔

For simplicity, we can take that the velocity vector at sections 1 and 2 is
normal to the surface, that is to say that in section 1: 𝑣⃗ = 𝑣1 (−𝑛̂1 ), and in section
2: 𝑣⃗ = 𝑣2 𝑛̂2 . Thus, Eq. 20 becomes:
𝑢 𝑢 𝑢
∫ 𝜌𝑔 ( + 𝐵) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 = − ∫ 𝜌𝑔 ( + 𝐵) 𝑣1 𝑑𝑆 + ∫ 𝜌𝑔 ( + 𝐵) 𝑣2 𝑑𝑆 (21)
𝑆𝐶 𝑔 𝑆1 𝑔 𝑆2 𝑔

:: 12 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑆2
𝑛̂2

𝑣⃗2
𝑆1

𝑆𝑚
𝑛̂1
SECTION 2

𝑣⃗1

SECTION 1

Fig. 9.- Application of the energy general equation to a streamtube

Now, we will do a strong assumption (but we will be able to overcome it,


so it will not be a problem in the future): let’s assume that the flow properties
are constant at each section of the stream flow (they are constant in the section,
but can change from one section to another). If this is the case, 𝑢, 𝑝⁄𝜌 and 𝑣, does
not depend on 𝑑𝑆 and we can take the 𝜌𝑔(𝑢 + 𝐵)𝑣 out of the integral:
𝑢 𝑢1 𝑢2
∫ 𝜌𝑔 ( + 𝐵) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 = −𝜌𝑔 ( + 𝐵1 ) 𝑣1 ∫ 𝑑𝑆 + 𝜌𝑔 ( + 𝐵2 ) 𝑣2 ∫ 𝑑𝑆 (22)
𝑆𝐶 𝑔 𝑔 𝑆1 𝑔 𝑆2

𝑢 𝑢1 𝑢2
∫ 𝜌𝑔 ( + 𝐵) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 = −𝜌𝑔 ( + 𝐵1 ) 𝑣1 𝑆1 + 𝜌𝑔 ( + 𝐵2 ) 𝑣2 𝑆2 (23)
𝑆𝐶 𝑔 𝑔 𝑔

But 𝑣1 𝑆1 is the discharge at section 1 which, by continuity, is equal to the


discharge at section 2, i.e.: 𝑣1 𝑆1 = 𝑣2 𝑆2 = 𝑄, and the integral over the surfaces is
reduced to a simpler expression:
𝑢 𝑢2 − 𝑢1
∫ 𝜌𝑔 ( + 𝐵) 𝑣⃗ ∙ 𝑛̂𝑑𝑆 = 𝜌𝑔𝑄 ( + 𝐵2 − 𝐵1 ) (24)
𝑆𝐶 𝑔 𝑔

Thus, for a steady flow through a stream tube with one entrance and one
exit (as in Fig. 9), the general equation of the energy is written as:
𝑢2 − 𝑢1 𝑑𝑄̂ 𝑑𝑊𝐸 𝑑𝑊𝜏
𝜌𝑔𝑄 ( + 𝐵2 − 𝐵1 ) = − − (25)
𝑔 𝑑𝑡 𝑑𝑡 𝑑𝑡

For flows of liquids (the most common in civil and environmental


engineering), the variation of internal energy is negligible, and we can write:

:: 13 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑑𝑄̂ 𝑑𝑊𝐸 𝑑𝑊𝜏


𝜌𝑔𝑄(𝐵2 − 𝐵1 ) = − − (26)
𝑑𝑡 𝑑𝑡 𝑑𝑡

For adiabatic flows (no heat exchange) of ideal fluids without work done
by the fluid, we obtain that the Bernoulli remains constant. :

𝐵2 = 𝐵1 = 𝐶𝑜𝑛𝑠𝑡. (27)

At this point, an objection could be made to the development of Eq. 26


when it is applied to real fluids. If the stream tube is considered to be a pipe, due
to the fluid viscosity, the velocity of the fluid particles in contact with the wall of
the pipe is zero and maximum at the axis of the pipe, violating the assumption
of constant velocity made to take the Bernoulli out of the surface integral and
obtain Eq. 25, which contains a term associated to the shear stress. In order to
avoid this inconsistency, the Bernoulli equation (Eq. 19) needs to be modified.
Let’s consider a flow with parallel streamlines (as the flow in a pipe of constant
diameter or the uniform flow in a channel) and analyse the term that contains 𝐵
in Eq. 18. As it was said, the variation of internal energy in liquids can be
neglected and we will leave it out. In a given section, we have:
𝑝 𝑣2
∫ 𝜌𝑔𝐵𝑣⃗ ∙ 𝑛̂𝑑𝑆 = ∫ 𝜌𝑔 (𝑧 + + ) 𝑣𝑑𝑆 (28)
𝑆 𝑆 𝜌𝑔 2𝑔

In any flow with parallel streamlines, the sum 𝑧 + 𝑝⁄(𝜌𝑔) remains


constant. In particular, we can evaluate that sum at the gravity centre (𝐺) of the
section:
𝑝 𝑣2 𝑝 𝑣2
∫ 𝜌𝑔 (𝑧 + + ) 𝑣𝑑𝑆 = ∫ 𝜌𝑔 (𝑧 + ) 𝑣𝑑𝑆 + ∫ 𝜌𝑔 𝑣𝑑𝑆 (29)
𝑆 𝜌𝑔 2𝑔 𝑆 𝜌𝑔 𝑆 2𝑔

𝑝 𝑣2 𝑝 𝑣2
∫ 𝜌𝑔 (𝑧 + + ) 𝑣𝑑𝑆 = 𝜌𝑔 (𝑧 + ) ∫ 𝑣𝑑𝑆 + 𝜌𝑔 ∫ 𝑣𝑑𝑆 (30)
𝑆 𝜌𝑔 2𝑔 𝜌𝑔 𝑆 𝑆 2𝑔

𝑝𝐺 𝑣2
∫ 𝜌𝑔𝐵𝑣𝑑𝑆 = 𝜌𝑔 (𝑧𝐺 + ) ∫ 𝑣𝑑𝑆 + 𝜌𝑔 ∫ 𝑣𝑑𝑆 (31)
𝑆 𝜌𝑔 𝑆 𝑆 2𝑔

As the velocity varies from zero to a maximum value, we want to express


the Bernoulli in terms of the mean velocity 𝑣̅ = 𝑄 ⁄𝐴, where 𝐴 is the area of the
surface 𝑆. Because 𝐵 depends on the square of 𝑣, the effect of using 𝑣̅ in the

:: 14 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Bernoulli is taking into account by means of a coefficient 𝛼 affecting the term


with𝑣̅ 2. Calling 𝐵 ∗ the Bernoulli based upon 𝑣̅ :
𝑝 𝑣̅ 2
𝐵∗ = 𝑧 + +𝛼 (32)
𝜌𝑔 2𝑔

We want to have:

∫ 𝜌𝑔𝐵𝑣𝑑𝑆 = ∫ 𝜌𝑔𝐵 ∗ 𝑣̅ 𝑑𝑆 (33)


𝑆 𝑆

Replacing Eq. 31 in the left side of Eq. 33 and Eq. 32 in the integral of the
right side:
𝑝𝐺 𝑣2 𝑝 𝑣̅ 2
𝜌𝑔 (𝑧𝐺 + ) ∫ 𝑣𝑑𝑆 + 𝜌𝑔 ∫ 𝑣𝑑𝑆 = ∫ 𝜌𝑔 (𝑧 + + 𝛼 ) 𝑣̅ 𝑑𝑆 (34)
𝜌𝑔 𝑆 𝑆 2𝑔 𝑆 𝜌𝑔 2𝑔

𝑝𝐺 𝑣2
𝜌𝑔 (𝑧𝐺 + ) ∫ 𝑣𝑑𝑆 + 𝜌𝑔 ∫ 𝑣𝑑𝑆
𝜌𝑔 𝑆 𝑆 2𝑔
𝑝𝐺 𝑣̅ 2 (35)
= 𝜌𝑔 (𝑧𝐺 + ) ∫ 𝑣̅ 𝑑𝑆 + 𝜌𝑔 ∫ 𝛼 𝑣̅ 𝑑𝑆
𝜌𝑔 𝑆 𝑆 2𝑔

But

∫ 𝑣𝑑𝑆 = ∫ 𝑣̅ 𝑑𝑆 = 𝑄 (36)
𝑆 𝑆

Using Eq. 36, Eq. 35 is reduced to:


𝑣2 𝑣̅ 2
∫ 𝑣𝑑𝑆 = ∫ 𝛼 𝑣̅ 𝑑𝑆 (37)
𝑆 2𝑔 𝑆 2𝑔

From where we can obtain an expression for the coefficient 𝛼 that takes
into account that the velocity is not constant in any section of the stream tube:

∫𝑆 𝑣 3 𝑑𝑆 (38)
𝛼=
𝑣̅ 3 𝐴
The coefficient 𝛼 is called Coriolis coefficient. It is easy to see that in order
to compute 𝛼 from Eq. 38, the velocity distribution has to be known, which is not
always possible, and many times it is determined from experiments or field
measurements. For turbulent flow in rectilinear pipes ~ 1.03 − 1.1 , depending
on 𝑅𝑒. As 𝛼 is close to 1, the value 𝛼 = 1 is commonly used. For a laminar flow
in a pipe, 𝛼 = 2.

:: 15 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Rigorously, Eq. 26 should be written as:

𝑑𝑄̂ 𝑑𝑊𝐸 𝑑𝑊𝜏


𝜌𝑔𝑄(𝐵2∗ − 𝐵1∗ ) = − − (39)
𝑑𝑡 𝑑𝑡 𝑑𝑡

with
𝑝 𝑣̅ 2

𝐵 =𝑧+ +𝛼 (32)
𝜌𝑔 2𝑔

However, it is customary that the bar ̅ over the velocity and the asterisk

of the Bernoulli are not written, and Eq. 32 becomes:
𝑝 𝑣2
𝐵=𝑧+ +𝛼 (40)
𝜌𝑔 2𝑔

As most of the applications deal with turbulent flow, 𝛼 = 1 is used, Eq. 32


ends written as:
𝑝 𝑣2
𝐵=𝑧+ + (41)
𝜌𝑔 2𝑔

A word of caution is necessary here: Eq. 41 is easily confused with Eq. 19,
although they are different. Eq. 41 can be applied to a real fluid with non-
uniform velocity profile, whereas Eq. 19 is restricted to uniform velocity profile,
something that cannot be attained by flows of real fluids with solid boundaries.
Following the common usage, we will write 𝐵, although we will be working
with𝐵 ∗ . Similarly, most of the time we will not write the Coriolis coefficient
because the flows will be turbulent.
If there is a hydraulic machine like a pump or a turbine, 𝑑𝑊𝐸 ⁄𝑑𝑡
corresponds to the power of the machine, 𝑃. We use to write the power divided
by (𝜌𝑔𝑄). For an ideal flow without heat exchange:
1 𝑑𝑊𝐸 𝑃
(𝐵2 − 𝐵1 ) = − =− (42)
𝜌𝑔𝑄 𝑑𝑡 𝜌𝑔𝑄

As we already mentioned, for a pump 𝑃 < 0 and for a turbine 𝑃 > 0.


Usually the notation ∆𝐵 = |𝑃|⁄(𝜌𝑔𝑄) is used, with the absolute value of the
power. The sign (+) for a turbine or (–) for a pump should be written explicitly in
the equation.
Let’s consider now the adiabatic flow of a real fluid, without external work.
The term − 𝑑𝑊𝜏 ⁄𝑑𝑡 is the power dissipated by viscous effects. It is customary in
hydraulics to work with the power divided by (𝜌𝑔𝑄). In this case, Eq. 39 becomes:

:: 16 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

1 𝑑𝑊𝜏
(𝐵2 − 𝐵1 ) = − (43)
𝜌𝑔𝑄 𝑑𝑡

The dissipated energy per unit weight of the fluid is denoted by, which
is defined as:
1 𝑑𝑊𝜏
≡ (44)
𝜌𝑔𝑄 𝑑𝑡

Thus, Eq. 43 is written as:

𝐵2 = 𝐵1 −  (45)

The energy loss is usually divided in two terms: one associated to friction
and other to singularities in the flow line. Let’s restrict our attention to the
friction loss. The energy loss per unit length is denoted by 𝐽, defined as:
𝑑𝐵
𝐽=− (46)
𝑑𝑥
Thus, the friction loss in a pipe of length 𝐿 is given by  = 𝐽𝐿.Computation
of 𝐽 is generally performed using the Darcy-Weisbach equation, which requires
to know the friction factor 𝑓. For a cylindrical pipe of diameter 𝐷:
1 𝑣2
𝐽=𝑓 (47)
𝐷 2𝑔

In general, 𝑓 depends of the flow regime, given by the Reynolds number,


𝑅𝑒, and the relative roughness of the pipe, 𝜀 ⁄𝐷, where 𝜀 is the size of the
roughness. It can be computed analytically for laminar flows, resulting 𝑓 =
64⁄𝑅𝑒. For turbulent flows, 𝑓 can be determined semi-analytically, as we will see
later. When Eq. 47 is applied to non-circular conduits, the diameter 𝐷 should be
replaced by the four times the hydraulic radius, defined as 𝑅𝐻 = 𝐴⁄ , where 𝐴
is the flow area and  is the wetted perimeter (i.e., 𝐷 = 4𝑅𝐻 )
Momentum theorem: It corresponds to the application of Newton’s second law
to a fluid system:
𝑑
(𝑚𝑣⃗) = 𝐹⃗ (48)
𝑑𝑡

For simplicity, let’s apply Eq. 48 to the same control volume used in the
derivation of the energy general equation (Fig. 9). It can be shown that the term

:: 17 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

that indicates the variation of momentum in Eq. 48, when apply to a fluid,
becomes:
𝑑 𝜕
(𝑚𝑣⃗) = ∫ 𝜌𝑣⃗ 𝑑𝑉 + ∫ 𝜌𝑣⃗ 𝑣⃗ ∙ 𝑛̂𝑑𝑆 (49)
𝑑𝑡 𝜕𝑡 𝑉𝐶 𝑆𝐶

The integral over the surfaces of the control volume is split in three terms:

∫ 𝜌𝑣⃗ 𝑣⃗ ∙ 𝑛̂𝑑𝑆 = ∫ 𝜌𝑣⃗ 𝑣⃗ ∙ 𝑛̂𝑑𝑆 + ∫ 𝜌𝑣⃗ 𝑣⃗ ∙ 𝑛̂𝑑𝑆 + ∫ 𝜌𝑣⃗ 𝑣⃗ ∙ 𝑛̂𝑑𝑆 (50)


𝑆𝐶 𝑆1 𝑆𝑚 𝑆2

Considering that there is not flow through the mantle and that 𝑣⃗1 ∙ 𝑛̂1 < 0,
𝑣⃗2 ∙ 𝑛̂2 > 0, constant properties at each transverse section of the stream tube, and
𝑣1 𝑆1 = 𝑣2 𝑆2 = 𝑄, Eq. 48 can be written as:
𝜕
∫ 𝜌𝑣⃗ 𝑑𝑉 + 𝜌𝑄(𝑣⃗2 − 𝑣⃗1 ) = 𝐹⃗ (51)
𝜕𝑡 𝑉𝐶

For a steady flow, Eq. 51 is simplified to:

𝜌𝑄(𝑣⃗2 − 𝑣⃗1 ) = 𝐹⃗ (52)

As it happened with the energy general equation, Eq. 52 is valid only for
uniform velocity profiles. In order to take into account the non-uniformity of the
velocity profile imposed by boundaries in flow of real fluids, the cross section
mean velocity is used in conjunction with a coefficient 𝛽 and Eq. 53 becomes:

𝜌𝑄(𝛽2 𝑣⃗2 − 𝛽1 𝑣⃗1 ) = 𝐹⃗ (52)

The coefficient 𝛽 is named Boussinesq coefficient and it is given by:

∫𝑆 𝑣 2 𝑑𝑆 (53)
𝛽=
𝑣̅ 2 𝐴
For most turbulent flows 𝛽 ~ 1.01 − 1.04, and the value 𝛽 = 1 is
considered. For a laminar flow in a pipe, 𝛽 = 4⁄3.
The uniform, steady, two-dimensional flow with free surface of a real
fluid: As an application of the three fundamental principles of the physics
applied to hydraulics, we will analyse the flow in an infinitely wide inclined
channel in which the permanent and uniform flow with a free surface exists. The
problem is sketched in Fig. 10. The bottom of the channel coincides with the 𝑥
axis which is inclined an angle 𝜃 with respect to the horizontal line. The flow

:: 18 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

depth is 𝐻 and the cross sectional mean velocity is 𝑈. There is not heat transfer
between the fluid sustem and the environment and the flow is turbulent.
The first thing that needs to be done when working with the equations
derived with the integral approach, is to choose the control volume. In this case
𝑉𝐶 corresponds to the volume defined by the segmented line. The application of
the continuity equation, energy equation and momentum theorem must be done
in this volume.

𝑝𝑎𝑡𝑚 𝑔⃗

𝑉𝐶
𝜃
Fig. 10.- Steady, uniform, 2-D free surface
flow
Continuity equation: As the flow is 2-D, we work with both the discharge and the
area of the flow section per unit width, i.e., 𝑞 = 𝑄 ⁄𝑏 and 𝑎 = 𝐴⁄𝑏. Calling (1) to
the entrance section and (2) to the exit section, and considering steady flow, the
application of Eq. 1 to the control volume reduces to 𝑞1 = 𝑞2 . Uniformity of the
flow means that the flow depth in both sections is the same, from where 𝑎1 =
𝑎2 = 𝐻. Thus, from Eq. 2: 𝑣1 = 𝑣2 = 𝑈.
Energy equation: The control volume corresponds to a steamtube, as shown in
Fig. 11. As the flow is turbulent we assume 𝛼 = 1. Application of Eq. 45:

𝐵2 = 𝐵1 −  (45)

(1)
(2)

𝑧1
𝑧2

Fig. 11.- Application of the energy equation to the control volume

:: 19 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

As the flow is uniform, the streamlines are parallels, and the sum 𝑧 +
𝑝⁄(𝜌𝑔) remains constant in any cross section. That means that the sum can be
evaluated at any 𝑧 of the section. For simplicity, we choose to evaluate it at the
free surface, because at that location the pressure is known (𝑝𝑎𝑡𝑚 ). Thus, Eq. 45
becomes:
𝑝𝑎𝑡𝑚 𝑈 2 𝑝𝑎𝑡𝑚 𝑈 2
𝑧2 + + = 𝑧1 + + − (54)
𝜌𝑔 2𝑔 𝜌𝑔 2𝑔

Thus:

 = 𝑧1 − 𝑧2 (55)

From the geometry:


𝑧1 − 𝑧2
sin 𝜃 = − (55)
∆𝑥

Using Eq. 46 and  = 𝐽∆𝑥, we obtain:

𝐽 = sin 𝜃 (56)

The result above obtained indicates that for steady, uniform open-channel
flows the gradient of the energy line is equal to the slope of the channel. Or,
equivalent, the energy line is parallel to the bottom of the channel (and parallel
to the free surface).
We can relate the flow velocity with the slope of the channel using Eq. 47.
It is easy to show that for this flow, 𝑅𝐻 = 𝐻. Thus, we have:
1 𝑈2
sin 𝜃 = 𝑓 (57)
4𝐻 2𝑔

The problem is not solved yet because we have not said anything about
the friction factor. We will get some relationships for it later.
Momentum theorem: In order to apply Newton’s second law to the control volume
of Fig. 10, we have to recognize that Eq. 52 is a vectorial equation, and it should
be applied to each component of the coordinate system. As we already have
defined the control volume, the first step in our analysis is to identify the forces
and momentum fluxes acting in 𝑉𝐶 . To do this, we will use Fig. 12.

:: 20 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

(1)
(2)
𝑔⃗
𝐹𝑝1
𝑈 𝐹𝑝2
𝑈
𝐹𝜏𝑜
𝑊

Fig. 12.- Application of the momentum theorem to the control volume

The 𝑦-component of the momentum equation is not relevant in this


analysis. It only says that the component of the weight in the 𝑦 direction is
equilibrated with the reaction normal to the bottom of the channel. Let’s analyse
the 𝑥-component of the momentum equation. Eq. 52 gives:

𝜌𝑄(𝛽2 𝑣2𝑥 − 𝛽1 𝑣1𝑥 ) = 𝐹𝑥 (58)

𝐹𝑥 is 𝑥-component the resultant of the forces acting in the control volume.


We can identify the following forces: 𝐹𝑝 , forces due to the fluid pressure; 𝑊,
weight of the fluid contained in the control volume; and 𝐹𝜏𝑜 , force due to the
friction between the fluid and the bottom. Thus:

𝐹𝑥 = 𝐹𝑝1 − 𝐹𝑝2 + 𝑊 sin 𝜃 − 𝐹𝜏0 (59)

As 𝑣1𝑥 = 𝑣2𝑥 = 𝑈, there is not variation of momentum, and Eq. 58 becomes


simply:

𝐹𝑝1 − 𝐹𝑝2 + 𝑊 sin 𝜃 − 𝐹𝜏0 = 0 (60)

We can evaluate the force due to the pressure as 𝐹𝑝 = 𝑝𝐺 𝐴, where 𝑝𝐺 is the


fluid pressure at the centre of gravity of the surface with area 𝐴. As the flow is
uniform, 𝐴 and 𝑝𝐺 are the same in sections (1) and (2), and the pressure does not
contribute in the equilibrium of forces given by Eq. 60. This indicates that for
the steady uniform flow in a channel, there is an equilibrium between the force
that generates the motion (the component of the weight in the flow direction)
and the force that opposes to the motion (due the friction at the wall):

:: 21 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑊 sin 𝜃 = 𝐹𝜏0 (61)

The weight (per unit width) is given by:

𝑊 = 𝜌𝑔𝐻∆𝑥 (62)

The force (per unit width) due to friction is:

𝐹𝜏0 = 𝜏0 ∆𝑥 (63)

In the above equation, 𝜏0 is the shear stress acting on the bottom of the
channel. Eqs. 61, 62 and 63 give:

𝜏0 = 𝜌𝑔𝐻 sin 𝜃 (64)

Using Eq. 56, we can write:

𝜏0 = 𝜌𝑔𝐻𝐽 (65)

Although Eq. 65 was obtained for steady and uniform flow, the same result
is obtained for steady and gradually varied flow.
Combining Eqs. 57 and 64:
1
𝜏0 = 𝜌𝑓𝑈 2 (66)
8

A “shear velocity”, 𝑢∗ , is defined as:


𝜏0
𝑢∗ = √ (67)
𝜌

From where:

𝑢∗ 𝑓
=√ (68)
𝑈 8

It is important to have in mind that 𝑢∗ is not a flow velocity. It is only a


way to write the friction stress acting on the wall, 𝜏0 . The term √𝜏0 ⁄𝜌 appears

:: 22 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

many times in the analysis of fluid flows. As it has dimensions of velocity, it is


called “friction velocity”.
Shear stress distribution in a steady, uniform, two-dimensional flow
with free surface of a real fluid: We obtained the shear stress acting on the
bottom of the channel (Eq. 64). As there is not shear applied on the free surface,
there the shear stress is zero. The question is: how does the shear stress changes
from zero at the free surface to 𝜏0 = 𝜌𝑔𝐻 sin 𝜃 on the bottom?. To answer this
question, we can applied the momentum theorem to the control volume of Fig.13.

𝑝𝑎𝑡𝑚 𝑔⃗

𝑉𝐶 𝜃
Fig. 13.- Control volume chosen to determine the variation of the shear
stress

The forces acting in 𝐶𝑉 are depicted in Fig. 14. Essentially, they are the
same than those shown in Fig. 12, but for a control volume that does not reach
the bottom of the channel.

(1)
(2)
𝐹𝑝𝑦1 𝑔⃗

𝑈 𝐹𝑝𝑦2
𝐹𝜏𝑦 𝑈
𝑊𝑦

Fig. 14.- Forces acting in the control volume

𝐹𝑝𝑦 is the force due to the fluid pressure acting on the surfaces comprised
between 𝐻 and 𝑦, 𝑊𝑦 is the weight of the fluid contained in the control volume,
and 𝐹𝜏𝑦 the force due to friction at a distance 𝑦 from the bottom. With the same
arguments that we obtained Eq. 61, we get:

:: 23 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑊𝑦 sin 𝜃 = 𝐹𝜏𝑦 (69)

The weight of fluid and force due to shear stress are given by:

𝑊𝑦 sin 𝜃 = 𝜌𝑔(𝐻 − 𝑦)∆𝑥 (70)

𝐹𝜏𝑦 = 𝜏𝑦 ∆𝑥 (71)

From where we obtain:

𝜏𝑦 = 𝜌𝑔(𝐻 − 𝑦) sin 𝜃 (72)

Eq. 72 indicates that the shear stress varies linearly with depth, from zero
on the free surface (𝑦 = 𝐻) to 𝜏𝑦 = 𝜌𝑔𝐻 sin 𝜃 on the bottom (𝑦 = 0), as sketched in
Fig. 15. Dividing Eq. 72 by Eq.64:
𝜏𝑦 𝑦
= (1 − ) (73)
𝜏0 𝐻

𝜏0 𝜏𝑦
Fig. 15.- Shear stress distribution

:: 24 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

BASIC EQUATIONS OF FLUID MECHANICS OBTAINED FROM THE


DIFFERENTIAL APPROACH

As we have seen, the equations obtained using the integral approach do


not provide information about the flow characteristics in a point of the flow
domain, but it considers mean velocities, and total forces acting on the control
volume or its surfaces. In order to take into account the variation of the velocity
in a given section, the coefficients 𝛼 and 𝛽 appears, but from their definition, the
velocity distribution is needed to compute them analytically. Now we will derive
the equations of continuity and momentum considering an infinitely small
control volume, which we will take to the limit that it becomes a point.

CONTINUITY EQUATION
Let’s analyse the flow of mass through an element of volume 𝑑𝑉 = 𝑑𝑥𝑑𝑦𝑑𝑧
immersed in the flow, as shown in Fig. 16. The velocity field is 𝑣⃗ = (𝑢, 𝑣, 𝑤).
Through this imaginary volume, the flow (represented by the streamlines in the
figure) passes transporting mass of fluid. Conservation of mass indicates that
the variation of mass per unit time inside of 𝑑𝑉 is equal to the mass rate that
enters to the volume, less that the mass rate that exit it. This statement can be

𝑑𝑧
𝑦
𝑑𝑦

𝑥 𝑑𝑥

Fig. 16.- Imaginary element of volume in the flow domain

written as:

𝜕𝑚
= 𝐺𝑖 − 𝐺𝑜 (74)
𝜕𝑡
where 𝑚 is the mass of fluid contained in 𝑑𝑉, given by

:: 25 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑑𝑚 = 𝜌𝑑𝑉 (75)

𝐺𝑖 is the mass rate that enters to 𝑑𝑉, and 𝐺𝑜 is the mass rate that exits the
volume. The fluid can enters (and exits) through any of the six surfaces that
define the volume. For the flow in the 𝑥 direction, we have a surface located at 𝑥
with normal (−𝑖̂), and other located at 𝑥 + 𝑑𝑥, with normal (+𝑖̂).The same
happens for the other coordinates. Thus, we can decompose the mass that enters
to 𝑑𝑉 in three terms, and write:

𝐺𝑖 = 𝐺𝑥 |𝑥 + 𝐺𝑦 |𝑦 + 𝐺𝑧 |𝑧 (76)

where the symbol |𝑥 stands for “evaluated at 𝑥”, and so on for the other
directions. Similarly, the mass that exits the volume can be written as:

𝐺𝑜 = 𝐺𝑥 |𝑥+𝑑𝑥 + 𝐺𝑦 |𝑦+𝑑𝑦 + 𝐺𝑧 |𝑧+𝑑𝑧 (77)

Replacing in Eq. 74:

𝜕
(𝜌𝑑𝑉) = 𝐺𝑥 |𝑥 + 𝐺𝑦 | + 𝐺𝑧 |𝑧 − (𝐺𝑥 |𝑥+𝑑𝑥 + 𝐺𝑦 | + 𝐺𝑧 |𝑧+𝑑𝑧 ) (78)
𝜕𝑡 𝑦 𝑦+𝑑𝑦

The terms evaluated at 𝑥 + 𝑑𝑥, 𝑦 + 𝑑𝑦 and in 𝑧 + 𝑑𝑧 can be related to those


evaluated at 𝑥, 𝑦 and in 𝑧 by means of a Taylor’s expansion. For example, the
term 𝐺𝑥 |𝑥+𝑑𝑥 is expanded as:
𝜕𝐺𝑥 1 𝜕 2 𝐺𝑥
𝐺𝑥 |𝑥+𝑑𝑥 = 𝐺𝑥 |𝑥 + | 𝑑𝑥 + | (𝑑𝑥)2 + ⋯ (79)
𝜕𝑥 𝑥 2 𝜕𝑥 2 𝑥

As 𝑑𝑥 is infinitely small, the term (𝑑𝑥)2 and higher powers of 𝑑𝑥 can be


neglected. Thus, the net flow of mass along the 𝑥 direction is:

𝜕𝐺𝑥
𝐺𝑥 |𝑥 − 𝐺𝑥 |𝑥+𝑑𝑥 = − | 𝑑𝑥 (80)
𝜕𝑥 𝑥

The mass rate is defined as the discharge times density. The discharge in
the 𝑥 direction is 𝑢𝑑𝐴𝑥 , where 𝑑𝐴𝑥 is the element of area of the surface with
normal 𝑖̂, 𝑑𝐴𝑥 = 𝑑𝑦𝑑𝑧. Thus, Eq. 80 becomes:

𝜕𝐺𝑥 𝜕
− | 𝑑𝑥 = − (𝜌𝑢)| 𝑑𝑥𝑑𝑦𝑑𝑧 (81)
𝜕𝑥 𝑥 𝜕𝑥 𝑥

In the same way:

:: 26 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕𝐺𝑦 𝜕 𝜕𝐺𝑧 𝜕
− | 𝑑𝑦 = − (𝜌𝑣)| 𝑑𝑥𝑑𝑦𝑑𝑧 − | 𝑑𝑧 = − (𝜌𝑤)| 𝑑𝑥𝑑𝑦𝑑𝑧 (82)
𝜕𝑦 𝑦 𝜕𝑦 𝑦
𝜕𝑧 𝑧 𝜕𝑧 𝑧

and Eq. 78 becomes:

𝜕𝜌 𝜕 𝜕 𝜕
𝑑𝑥𝑑𝑦𝑑𝑧 = − ( (𝜌𝑢) + (𝜌𝑣) + (𝜌𝑤)) 𝑑𝑥𝑑𝑦𝑑𝑧 (83)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧

Thus, the continuity equation in the differential approach is written as:

𝜕𝜌 𝜕 𝜕 𝜕
+ (𝜌𝑢) + (𝜌𝑣) + (𝜌𝑤) = 0 (84)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧

or, in vectorial form:

𝜕𝜌
+ ∇ ∙ (𝜌𝑣⃗) = 0 (85)
𝜕𝑡
For an incompressible fluid, Eq. 85 is greatly reduced to:

∇ ∙ 𝑣⃗ = 0 (86)

Or, equivalently:

𝜕𝑢 𝜕𝑣 𝜕𝑤
+ + =0 (87)
𝜕𝑥 𝜕𝑦 𝜕𝑧

Eq. 87 (or 86) is the equivalent to Eq.1 of the integral approach.

MOMENTUM EQUATION

Before applying the Newton’s second law to the elementary volume of Fig.
16, it is important to remember that two kind of forces can be identified in the
element of fluid: Forces that depend on the amount of matter and that their
application point is the centre of gravity of the element of fluid (for example,
weight). They are called “body forces”, and notated as 𝐹⃗𝐵 . The second kind of
forces act on the surfaces that define the volume of fluid (like pressure and shear
stresses), and “surface forces”, and written as 𝐹⃗𝑆 . Thus, Newton’s second law can
be written as:

:: 27 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑑
(𝑚𝑣⃗) = 𝐹⃗𝑆 + 𝐹⃗𝐵 (88)
𝑑𝑡

As the total mass of the system remains constant, Eq. 88 evolves into:
𝑑𝑣⃗
𝜌𝑑𝑉 = 𝐹⃗𝑆 + 𝐹⃗𝐵 (89)
𝑑𝑡

The total or material derivative of the velocity can be expressed as:


𝑑𝑣⃗ 𝜕𝑣⃗ 𝜕𝑣⃗ 𝜕𝑣⃗ 𝜕𝑣⃗ 𝜕𝑣⃗
= + (𝑣⃗ ∙ ∇)𝑣⃗ = +𝑢 +𝑣 +𝑤 (91)
𝑑𝑡 𝜕𝑡 𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧

The body force depends of the amount of matter, i.e., it is proportional to


the mass of fluid in the 𝑑𝑉. As the mass is given by 𝜌𝑑𝑉, we can express the body
force in terms of a body per unit mass, 𝑓⃗𝐵 :

𝐹⃗𝐵 = 𝜌𝑓⃗𝐵 𝑑𝑉 (92)

The most common body force is that due to gravity. In this case, 𝑓⃗𝐵 = 𝑔⃗,
and Eq. 92 becomes:

𝐹⃗𝐵 = 𝜌𝑔⃗𝑑𝑉 (93)

The analysis of surface forces requires to evaluate the forces acting on the
six surfaces that define de element of fluid volume, as shown in Fig. 17. There
are three surface forces acting on each surface. The surface stresses are
represented as 𝜏𝑖𝑗 , where 𝑖 indicates the direction of the vector normal to the
surface and 𝑗 the direction of the force (where 𝑖 and 𝑗 can take the values 𝑖̂, 𝑗̂ or
𝑘̂ ). Thus, the net force in the direction 𝑗 is the result of the stresses of the
surfaces located at 𝑥, with normal (−𝑖̂); at 𝑥 + 𝑑𝑥, with normal (+𝑖̂); at 𝑦, with
normal (−𝑗̂); at 𝑦 + 𝑑𝑦, with normal (+𝑗̂); and at 𝑧, with normal (−𝑘̂); at 𝑧 + 𝑑𝑧,
with normal (+𝑘̂). Obviously, the direction 𝑗 can be 𝑖̂, 𝑗̂ or 𝑘̂.
The net surface force has components in the three directions:

𝐹⃗𝑆 = 𝐹𝑆𝑥 𝑖̂ + 𝐹𝑆𝑦 𝑗̂ + 𝐹𝑆𝑧 𝑘̂ (94)

As it was said, six surfaces contribute to each component of Eq. 94. For
the component in the 𝑥 direction we have:

:: 28 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝐹𝑆𝑥 = 𝐹𝑆𝑥𝑥 |𝑥 + 𝐹𝑆𝑥𝑥 |𝑥+𝑑𝑥 + 𝐹𝑆𝑦𝑥 |𝑦 + 𝐹𝑆𝑦𝑥 |𝑦+𝑑𝑦 + 𝐹𝑆𝑧𝑥 |𝑧 + 𝐹𝑆𝑧𝑥 |𝑧+𝑑𝑧 (95)

In term of the stresses, any surface force can be written as 𝐹𝑆𝑖𝑗 = 𝜏𝑖𝑗 𝑑𝑆⃗𝑗 ,
where 𝑑𝑆⃗𝑗 is the surface on which the stress is acting, with a vector normal to
the surface in the 𝑗 direction. Thus, Eq. 95 is given by:

𝜏𝑧𝑧
𝜏𝑧𝑦

𝜏𝑧𝑥 𝜏𝑥𝑧

𝜏𝑦𝑧 𝜏𝑥𝑦

𝑧 𝜏𝑥𝑥
𝑑𝑧
𝜏𝑦𝑦 𝜏𝑦𝑥
𝑦
𝑑𝑦

𝑥 𝑑𝑥

Fig. 17.- Surface forces acting on the three visible surfaces of


the element of fluid volume

𝐹𝑆𝑥 = 𝜏𝑥𝑥 |𝑥 (−𝑖̂)𝑑𝑆𝑥 + 𝜏𝑥𝑥 |𝑥+𝑑𝑥 (𝑖̂)𝑑𝑆𝑥 + 𝜏𝑦𝑥 |𝑦 (−𝑗̂)𝑑𝑆𝑦 + 𝜏𝑦𝑥 |𝑦+𝑑𝑦 (𝑗̂)𝑑𝑦
(96)
+ 𝜏𝑧𝑥 |𝑧 (−𝑘̂)𝑑𝑆𝑧 + 𝜏𝑧𝑥 |𝑧+𝑑𝑧 (𝑘̂)𝑑𝑆𝑧

The surfaces on which the stresses are acting correspond to:

𝑑𝑆𝑥 = 𝑑𝑦𝑑𝑧 𝑑𝑆𝑦 = 𝑑𝑥𝑑𝑧 𝑑𝑆𝑧 = 𝑑𝑥𝑑𝑦 (97)

Expanding in Taylor’s series and neglecting the terms of second order and
higher, we have:

:: 29 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕𝜏𝑥𝑥
𝜏𝑥𝑥 |𝑥+𝑑𝑥 = 𝜏𝑥𝑥 |𝑥 + | 𝑑𝑥 + ⋯
𝜕𝑥 𝑥
𝜕𝜏𝑦𝑥
𝜏𝑦𝑥 |𝑦+𝑑𝑦 = 𝜏𝑦𝑥 |𝑦 + | 𝑑𝑦 + ⋯ (98)
𝜕𝑦 𝑦
𝜕𝜏𝑧𝑥
𝜏𝑧𝑥 |𝑧+𝑑𝑧 = 𝜏𝑧𝑥 |𝑧 + | 𝑑𝑧 + ⋯
𝜕𝑧 𝑧

Replacing Eqs. 97 and 98 in Eq. 96, we get:


𝜕𝜏𝑥𝑥 𝜕𝜏𝑦𝑥 𝜕𝜏𝑧𝑥
𝐹𝑆𝑥 = ( + + ) 𝑑𝑥𝑑𝑦𝑑𝑧 (99)
𝜕𝑥 𝜕𝑦 𝜕𝑧

Following the same procedure for the components along 𝑦 and 𝑧 of 𝐹⃗𝑆 , we
obtain:
𝜕𝜏𝑥𝑦 𝜕𝜏𝑦𝑦 𝜕𝜏𝑧𝑦
𝐹𝑆𝑦 = ( + + ) 𝑑𝑥𝑑𝑦𝑑𝑧 (100)
𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕𝜏𝑥𝑧 𝜕𝜏𝑦𝑧 𝜕𝜏𝑧𝑧


𝐹𝑆𝑧 = ( + + ) 𝑑𝑥𝑑𝑦𝑑𝑧 (101)
𝜕𝑥 𝜕𝑦 𝜕𝑧

Replacing Eqs. 91, 93, 99, 100 and 101 in Eq. 94, Newton’s second law is
written as:

𝜕𝑣⃗ 𝜕𝑣⃗ 𝜕𝑣⃗ 𝜕𝑣⃗


𝜌( +𝑢 +𝑣 + 𝑤 ) = 𝜌𝑔⃗ +
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧
(102)
𝜕𝜏𝑥𝑥 𝜕𝜏𝑦𝑥 𝜕𝜏𝑧𝑥 𝜕𝜏𝑥𝑦 𝜕𝜏𝑦𝑦 𝜕𝜏𝑧𝑦 𝜕𝜏𝑥𝑧 𝜕𝜏𝑦𝑧 𝜕𝜏𝑧𝑧
( + + ) 𝑖̂ + ( + + ) 𝑗̂ + ( + + ) 𝑘̂
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

We can recognize that 𝜏𝑥𝑥 , 𝜏𝑦𝑥 , 𝜏𝑧𝑥 , 𝜏𝑥𝑦 , … are the elements of matrix array.
Actually, they are a tensor, and 𝜏𝑖𝑗 corresponds to the stress tensor:

𝜏𝑥𝑥 𝜏𝑦𝑥 𝜏𝑧𝑥


𝜏𝑖𝑗 = (𝜏𝑥𝑦 𝜏𝑦𝑦 𝜏𝑧𝑦 ) (103)
𝜏𝑥𝑧 𝜏𝑦𝑧 𝜏𝑧𝑧

Eq. 102 is greatly simplified when it is written in vectorial form:


𝜕𝑣⃗
𝜌 ( + (𝑣⃗ ∙ ∇)𝑣⃗) = 𝜌𝑔⃗ + ∇ ∙ 𝜏𝑖𝑗 (104)
𝜕𝑡

:: 30 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

As Eq. 104 is vectorial, we can write three equations, one for each
component.
Component 𝑥:
𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝜏𝑥𝑥 𝜕𝜏𝑦𝑥 𝜕𝜏𝑧𝑥
𝜌( +𝑢 +𝑣 + 𝑤 ) = 𝜌𝑔𝑥 + + +
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧 (105)

Component 𝑦:
𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝜏𝑥𝑦 𝜕𝜏𝑦𝑦 𝜕𝜏𝑧𝑦
𝜌( +𝑢 +𝑣 + 𝑤 ) = 𝜌𝑔𝑦 + + +
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧 (106)

Component 𝑧:
𝜕𝑤 𝜕𝑤 𝜕𝑤 𝜕𝑤 𝜕𝜏𝑥𝑧 𝜕𝜏𝑦𝑧 𝜕𝜏𝑧𝑧
𝜌( +𝑢 +𝑣 +𝑤 ) = 𝜌𝑔𝑧 + + +
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧 (107)

In the above equations, 𝑔⃗ = (𝑔𝑥 , 𝑔𝑦 , 𝑔𝑧 ) was used.


Eq. 104 (or Eqs. 105 to 107) are the Cauchy’s equations. It has to be noted
that up to now, we have not imposed explicitly at what kind of matter is applied
(the only restriction is that it has to be continuous, in a gravitational field). It
can be applied as much as solids as to fluids. We define the kind of matter
through the so-called “constitutive relationships”. They are relationships
between the stresses 𝜏𝑖𝑗 and the deformation or deformation rate of the matter.
We will restrict the analysis to Newtonian fluids (for example fluids like water,
air, oil, etc.). We will not perform the derivation of the constitutive relationships
for Newtonian fluids but give the final result, obtained by Stokes in 1845.
Basically, the assumptions made by Stokes to derive the equation of motion of
Newtonian fluids are the following:
1. The fluid is a continuum
2. If there is not motion, the equations of hydrostatic should be recovered
3. There is, at the most, a linear relationship between the stresses and
the angular deformation rate of an element of fluid
4. The fluid is isotropic, i.e., the constitutive relationship is independent
of the direction or coordinate system.
The second assumption indicates that if there is no motion, the angular
deformation rates should be zero and the normal stresses acting on the element
of fluid should reduce to the pressure. This means that

𝜏𝑖𝑗 = −𝑝𝛿𝑖𝑗 if 𝑣⃗ = 0 (108)

:: 31 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝛿𝑖𝑗 is the Kronecker’s delta, defined as:


1 if 𝑖 = 𝑗
𝛿𝑖𝑗 = { (109)
0 if 𝑖 ≠ 𝑗
The sign (−) in Eq. 108 results from the fact that the pressure points
towards de element of fluid. The final result from Stokes analysis is:
𝜕𝑢𝑖 𝜕𝑢𝑗
𝜏𝑖𝑗 = −𝑝𝛿𝑖𝑗 + 𝜇 ( + ) + 𝛿𝑖𝑗 ∇ ∙ 𝑣⃗ (110)
𝜕𝑥𝑗 𝜕𝑥𝑖

(As before, 𝑖 and 𝑗 indicate the component along 𝑥, 𝑦 or 𝑧)


𝜇 is the dynamic viscosity and  is the second coefficient of viscosity.
Stokes defined the mechanical pressure, 𝑝̅ , as the negative average of the
normal stresses:
𝜏𝑥𝑥 + 𝜏𝑦𝑦 + 𝜏𝑧𝑧
𝑝̅ = − (111)
3

From Eq. 110, we get


𝜕𝑢
𝜏𝑥𝑥 = −𝑝 + 𝜇2 + ∇ ∙ 𝑣⃗
𝜕𝑥
𝜕𝑣
𝜏𝑦𝑦 = −𝑝 + 𝜇2 + ∇ ∙ 𝑣⃗ (112)
𝜕𝑦
𝜕𝑤
𝜏𝑧𝑧 = −𝑝 + 𝜇2 + ∇ ∙ 𝑣⃗
𝜕𝑧
Replacing the last equations in Eq. 111:
2
𝑝̅ = 𝑝 − ( 𝜇 + ) ∇ ∙ 𝑣⃗ (113)
3

Note that Eq. 113 gives an interesting result: the mechanical pressure, 𝑝̅ ,
is not equal to the thermodynamic pressure, 𝑝. An interesting issue is raised
with the second coefficient of viscosity. In his work, Stokes assumed that
(2𝜇 ⁄3 + ) = 0, meaning that  < 0. However, some measurements indicate that
 > 0. Nevertheless, the value of  should not bother us because for
incompressible fluids ∇ ∙ 𝑣⃗ = 0 (Eq. 86), and in this case 𝑝̅ = 𝑝. Also, for
uncompressible fluids, Eq. 110 is reduced to:
𝜕𝑢𝑖 𝜕𝑢𝑗
𝜏𝑖𝑗 = −𝑝𝛿𝑖𝑗 + 𝜇 ( + ) (114)
𝜕𝑥𝑗 𝜕𝑥𝑖

Thus, the nine components of the tensor 𝜏𝑖𝑗 (Eq. 103) are:

:: 32 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕𝑢
𝜏𝑥𝑥 = −𝑝 + 𝜇2 (115)
𝜕𝑥

𝜕𝑣
𝜏𝑦𝑦 = −𝑝 + 𝜇2 (116)
𝜕𝑦

𝜕𝑤
𝜏𝑧𝑧 = −𝑝 + 𝜇2 (117)
𝜕𝑧

𝜕𝑢 𝜕𝑣
𝜏𝑥𝑦 = 𝜏𝑦𝑥 = 𝜇 ( + ) (118)
𝜕𝑦 𝜕𝑥

𝜕𝑢 𝜕𝑤
𝜏𝑥𝑧 = 𝜏𝑧𝑥 = 𝜇 ( + ) (119)
𝜕𝑧 𝜕𝑥

𝜕𝑣 𝜕𝑤
𝜏𝑦𝑧 = 𝜏𝑧𝑦 = 𝜇 ( + ) (120)
𝜕𝑧 𝜕𝑦

The terms containing 𝜇 are the viscous stresses. Introducing Eqs. 115 to 120 in
Eqs.105 to 107 and using Eq. 87, we obtain the momentum equations for an
incompressible Newtonian fluid with constant viscosity in the gravitational field:
Component 𝑥:
𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑝 𝜕 2𝑢 𝜕 2𝑢 𝜕 2𝑢
𝜌( +𝑢 +𝑣 +𝑤 )=− + 𝜇 ( 2 + 2 + 2 ) + 𝜌𝑔𝑥 (121)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧

Component 𝑦:
𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑝 𝜕 2𝑣 𝜕 2𝑣 𝜕 2𝑣
𝜌( +𝑢 +𝑣 +𝑤 )=− + 𝜇 ( 2 + 2 + 2 ) + 𝜌𝑔𝑦 (122)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑧

Component 𝑧:
𝜕𝑤 𝜕𝑤 𝜕𝑤 𝜕𝑤 𝜕𝑝 𝜕 2𝑤 𝜕 2𝑤 𝜕 2𝑤
𝜌( +𝑢 +𝑣 +𝑤 )=− +𝜇( 2 + + ) + 𝜌𝑔𝑧 (123)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑥 𝜕𝑦 2 𝜕𝑧 2

Eqs. 121, 122 and 123 can be written in a more compact form using
vectorial notation:
𝜕𝑣⃗
𝜌 ( + (𝑣⃗ ∙ ∇)𝑣⃗) = −∇𝑝 + 𝜇∇2 𝑣⃗ + 𝜌𝑔⃗ (124)
𝜕𝑡

Eqs 121 to 123 (or Eq. 124) are called the Navier-Stokes equations. The
steps given in this note to obtain the equations are based in the work of George

:: 33 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Gabriel Stokes “On the Theories of Internal Friction of Fluids in Motion”,


published in 1845, but in 1822 Claude-Louis Marie Henri Navier presented the
article “Memoire sur les lois du mouvement des fluides” (published in 1823),
where he got the same sets of Eqs. 121 to 123, but instead of the viscosity 𝜇
multiplying the second derivatives of the velocity, he had a coefficient −𝜀 that
arose from the resistance generated by the slip of adjacent layers of fluid
molecules (Navier, 1823, p. 416). Although the equations obtained by Navier
have the right form, he was not able to link the origin of the molecular forces
with viscosity.
The unknowns of a problem of fluid motion are four: the three components
of the velocity field (𝑢, 𝑣, 𝑤) and the pressure (𝑝). In general, they are functions
of the space and time: (𝑥, 𝑦, 𝑧, 𝑡).Conceptually, the problem is already solved,
because we have a system of four partial differential equations, with their
corresponding boundary and initial conditions. Three equations corresponds to
the Navier-Stokes equations (Eqs. 121, 122 and 123) and the fourth is the
continuity equation (Eq. 87). However, the set differential equations is highly
complex, due to the nonlinear terms (𝑢 𝜕𝑢⁄𝜕𝑥 , 𝑣 𝜕𝑢⁄𝜕𝑦 , … ) of the momentum
equations. Thus, analytical solutions are restricted to a few simple cases. The
existence of solutions of the Navier-Stokes equations and if they are unique (for
the general case) has not been proved yet, and it is one of the problem of the
millennium. There is a prize of US$ 1 million to be awarded to whom can prove
that the solution exists and it is unique. More information about the prize can
be found in the website www.claymath.org/millennium-problems/millennium-
prize-problems. As an anecdote, we can mention that in 2013, Mukhtarbay
Otelbaev, from the Institute of Mathematics and Mathematical Modelling,
Kazakhastan, published the paper (in Russian) “Existence of a strong solution
of the Navier-Stokes equation” in the Mathematical Journal (Almaty), Vol 13,
No. 4, pp. 5-104, where through 101 pages he would have solved the problem.
However, about a month later, Otelbaev recognized that he made a mistake in
one of the steps of his deduction (Moskvitch, 2014). As a recognition of the
importance of the Navier-Stokes, one student of the Department of Mining
Engineering of the University of Chile tattooed the equations in his arm, as
shown in Fig. 18xx. As he works with slurries that behave as non-Newtonian
fluid, the constitutive relationship differs from Eqs. 115 to 120 and he wrote only
𝜏. For non-Newtonian fluids, the expression for 𝜏 can be very complex and the
student will need his other arm to tattoo them!.
In the development of the Navier-Stokes equations we have not imposed
any restriction regarding the flow regime, and they are valid as much for laminar
as for turbulent flows. However, analytical solutions are obtained only for
laminar flows. The approach to turbulent flows by means of the Navier-Stokes

:: 34 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

equations is done numerically, and it has open a complete field of research y fluid
mechanics, the computational fluid dynamics (CFD).

Fig. 18xx.- Tattoo in the arm of one student of the Mining Engineering
Department of the University of Chile.

The steady, uniform 2-D laminar flow with free surface over an inclined
plane. To fix ideas, we will solve now the problem that we already analysed
using the integral approach. A fluid of density 𝜌 and dynamic viscosity 𝜇 flows
over a plane inclined an angle 𝜃 with a flow depth 𝐻 due to the action of gravity.
The flow is steady, uniform and laminar. The problem is to determine the
velocity and pressure distributions.
We choose the coordinate system indicated in Fig. 18. As the flow is 2-D,
we can omit the 𝑧 component of the momentum equation and drop 𝑤 and 𝑧

𝑝𝑎𝑡𝑚 𝑔⃗

𝜃
Fig. 18.- Steady, uniform, 2-D free surface laminar flow over an inclined
plane
derivatives. Thus the equation of continuity and Navier-Stokes are reduced to:

:: 35 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Continuity equation:

𝜕𝑢 𝜕𝑣
+ =0 (125)
𝜕𝑥 𝜕𝑦

As the flow is uniform, the streamlines are parallel to the 𝑥 axis (the free
surface is parallel to the bottom). Thus, there is not component of the velocity in
the 𝑦 direction. Therefore:

𝑣=0 0≤𝑦≤𝐻 (126)

and

𝜕𝑢
=0 (127)
𝜕𝑥
Navier-Stokes equation in the 𝑥 direction:
𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑝 𝜕 2𝑢 𝜕 2𝑢
𝜌( +𝑢 +𝑣 )=− + 𝜇 ( 2 + 2 ) + 𝜌𝑔𝑥 (128)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦

Navier-Stokes equation in the 𝑦 direction:


𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑝 𝜕 2𝑣 𝜕 2𝑣
𝜌( +𝑢 +𝑣 )=− + 𝜇 ( 2 + 2 ) + 𝜌𝑔𝑦 (129)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦

Eqs. 128 and 129 are greatly simplified with the conditions of the problem
and the result obtained from the continuity equation (Eqs. 126 and 127)
The condition of steady state means that the partial derivatives with
respect to 𝑡 are zero. This condition and Eqs. 126 and 127 leads to:
Navier-Stokes equation in the 𝑥 direction:
𝜕𝑝 𝜕 2𝑢
0=− + 𝜌𝑔𝑥 + 𝜇 2 (130)
𝜕𝑥 𝜕𝑦

Navier-Stokes equation in the 𝑦 direction:


𝜕𝑝
0=− + 𝜌𝑔𝑦 (131)
𝜕𝑦

The last equation indicates that the pressure distribution varies linearly
with 𝑦, i.e., we have a hydrostatic pressure distribution. Integrating Eq. 131:

:: 36 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑝 = 𝜌𝑔𝑦 𝑦 + 𝐶1 (𝑥) (132)

𝜕(Eq. 132)⁄𝜕𝑥 and replacing it in Eq. 130:


𝜕𝐶1 𝜕 2𝑢
0=− + 𝜌𝑔𝑥 + 𝜇 2 (133)
𝜕𝑥 𝜕𝑦

From Eq.127 we know that the third term of Eq. 133 does not depends on
𝑥. It means that we can take that term as a constant 𝐾 if we integrate Eq. 133
with respect to 𝑥:
𝜕𝐶1
= 𝜌𝑔𝑥 + 𝐾 (134)
𝜕𝑥

𝐶1 = (𝜌𝑔𝑥 + 𝐾)𝑥 + 𝐶 (135)

𝐶 is a constant pure because we already know that 𝐶1 does not depend on 𝑦.


Replacing 𝐶1 in Eq. 132:

𝑝 = 𝜌𝑔𝑦 𝑦 + (𝜌𝑔𝑥 + 𝐾)𝑥 + 𝐶 (136)

To determine 𝐶 we have to apply the boundary condition for the pressure.


At the free surface we have atmospheric pressure. Working with relative
pressures, 𝑝𝑎𝑡𝑚 = 0 and the boundary condition becomes:

For any 𝑥 at 𝑦 = 𝐻, 𝑝 = 0 (137)

Thus:

𝐶 = −𝜌𝑔𝑦 𝐻 − (𝜌𝑔𝑥 + 𝐾)𝑥 (138)

Replacing now the 𝐶 in Eq. 136:

𝑝 = 𝜌𝑔𝑦 (𝑦 − 𝐻) (139)

From Fig. 18 it is easy to see that

𝑔𝑥 = 𝑔 sin 𝜃 , 𝑔𝑦 = −𝑔 cos 𝜃 (140)

:: 37 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Finally, the pressure distribution is:

𝑝 = 𝜌𝑔 cos 𝜃 (𝐻 − 𝑦) (141)

Known the pressure, we can go back to Eq. 130. From Eq. 141, we have
𝜕𝑝⁄𝜕𝑥 = 0, thus:
𝜕 2𝑢
𝜇 2 = −𝜌𝑔𝑥 (142)
𝜕𝑦

Using the cinematic viscosity,  = 𝜇 ⁄𝜌 and Eq. 140:

𝜕 2𝑢 𝑔
= − sin 𝜃 (143)
𝜕𝑦 2 

Integrating Eq. 143 twice with respect to 𝑦:

𝑔𝑦2
𝑢 = − sin 𝜃 + 𝐴𝑦 + 𝐵 (144)
 2

The constants 𝐴 and𝐵 are determined from the boundary conditions. At


the bottom we have the non-slip condition, i.e., the fluid velocity is the same that
the bottom velocity. It means:

For any 𝑥 at 𝑦 = 0, 𝑢=0 (145)

The second boundary condition is on the free surface. As there is not shear
stress applied on it, the condition is:

For any 𝑥 at 𝑦 = 𝐻, 𝜏𝑦𝑥 = 0 (146)

That is to say:
𝜏𝑦𝑥 | =𝜇
𝜕𝑢
| =0 (147)
𝑦=𝐻 𝜕𝑦 𝑦=𝐻

The following values are obtained:


𝑔
𝐴 = sin 𝜃 𝐻

(148)
𝐵=0

:: 38 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Finally, the velocity distribution is given by:


𝑔 𝑦2
𝑢= sin 𝜃 (𝐻𝑦 − ) (149)
 2

From Eq. 149 we see that the velocity distribution is parabolic, with a
maximum value at the free surface equal to 𝑢 = 𝑔 sin 𝜃 𝐻 2 ⁄(2). We can also
compute the average velocity, 𝑢̅ = 𝑈:
1 𝐻
𝑈 = ∫ 𝑢 𝑑𝑦
𝐻 0
(150)
1𝑔
𝑈= sin 𝜃 𝐻 2
3

We can write Eq. 150 in dimensionless form multiplying both sides by 𝑈⁄(𝑔𝐻):
𝑈 2 1 𝑈𝐻
= sin 𝜃 (151)
𝑔𝐻 3 

We recognize in the left hand side of Eq. 151 a Reynolds number based on the
flow depth:
𝑈𝐻
𝑅𝑒𝐻 =
 (152)

In the left hand side there is a dimensionless number that appears when we are
work when gravity forces are important. It is the square of Froude number, which is
defined as:
𝑈
𝐹𝑟 =
√𝑔𝐻 (153)

The equation that relates the average velocity of the flow with its depth is called
in hydraulics “resistance law”, which in this case can be written as:
1
𝐹𝑟 2 = 𝑅𝑒𝐻 sin 𝜃
3 (154)

As we know the velocity distribution, we can also compute the coefficients 𝛼 and
𝛽 introduced in the integral approach (Eqs. 38 and 53), resulting 𝛼 = 54⁄35 ≈ 1.543
and 𝛽 = 18⁄15 = 1.2.
We can also compute the shear stress distribution and its value at the bottom,
that we called 𝜏0 :

:: 39 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕𝑢
𝜏𝑦𝑥 = 𝜇
𝜕𝑦
(155)
𝜏𝑦𝑥 = 𝜌𝑔 sin 𝜃(𝐻 − 𝑦)

At the bottom (𝑦 = 0):

(156)
𝜏0 = 𝜌𝑔𝐻 sin 𝜃

Note that the relationships given by Eqs. 155 and 156 are the same than
those obtained with the application of the momentum theorem (Eqs. 72 and 64).
Something that we cannot obtain from the integral approach is the resistance
relationship given by Eq. 154. The equivalent equation in the integral approach
is the Darcy-Weisbach equation (Eq. 47 that evolved into Eq. 57). However, this
equation requires to know the friction factor 𝑓, that should be determined from
other way (theoretically, numerically or experimentally). Using Eqs. 57 and 154
we can obtain the friction factor for a steady, uniform, 2-D laminar flow in a
channel:
𝑈2 𝐹𝑟 2
sin 𝜃 = 𝑓 =3 (157)
8𝑔𝐻 𝑅𝑒𝐻

From where:
24
𝑓=
𝑅𝑒𝐻 (158)

Do not mistake the expression given by Eq. 158 with the expression for a
cylindrical pipe, 𝑓 = 64⁄𝑅𝑒 (and shown graphically in Moody’s diagram).
Replacing in 𝑅𝑒 the diameter 𝐷 by 4𝑅𝐻 we do not obtain Eq. 158 but 𝑓 = 16⁄𝑅𝑒𝐻 ,
which is wrong. The factor 24 has been confirmed experimentally (Chow, 1988).

:: 40 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

2.- REYNOLDS’ EQUATIONS FOR THE TRBULENCE

Turbulent flows are unsteady flows, in the sense that at any point of the
flow domain, the flow properties (𝑣⃗, 𝑝) are fluctuating in time. This is shown in
Fig.2.1, where the three components of the velocity measured with an acoustic
Dopper velocimeter (ADV) are presented. It is easily seen the fluctuating
characteristic of the data. It is rather difficult to use the data as it is presented
in the figure. Osborne Reynolds in 1895 proposed that in the turbulent flow
velocity and pressure could be decomposed in two terms, called by him mean-
mean-motion and relative-mean-motion (Reynolds, 1895). The peer reviewers of
the paper were G.G. Stokes and H. Lamb. In the first response that they sent to
the editor (Lord Rayleigh), Stokes acknowledged that he did not understand the
work, and Lamb indicated that much of the paper was obscure (Jackson and
Launder, 2007). Currently, we call those terms mean (or average) and
fluctuation. But, mean or average of what?

VELOCITY TIME SERIES (z = 17.33 cm)


4
u(t)
3.5

2.5
(cm/s)

1.5
v(t)
1

0.5

0
w(t)
-0.5
0 20 40 60 80 100 120 140 160

t (seg)
Fig. 2.1.- Record of the three components of the velocity measured at
one location in a turbulent flow in an open channel

:: 41 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

The simplest (or easier to understand) is the temporal mean or temporal


average. Considering, for example, the record along the time 𝑡 of the 𝑢 component
of the velocity at a given location 𝑥⃗ of the flow domain, the temporal average is
defined as:

1 𝑡0 +𝑇
𝑢̅(𝑥⃗) = ∫ 𝑢(𝑥⃗, 𝑡) 𝑑𝑡 (2.1)
𝑇 𝑡0

WARNING! Unfortunately, we are using the same notation for to different


averages. We used the overbar ̅ previously to denote the average velocity in a
cross section (𝑣̅ = 𝑄/𝐴), and now we are using the overbar to denote an average
on time. Be careful!
If the average given by Eq. 2.1 does not depend on 𝑡0 , the process is called
statistically stationary.
According to the Reynolds decomposition, the velocity 𝑢(𝑥⃗, 𝑡) can be split
in two components: its temporal average 𝑢̅(𝑥⃗) and a fluctuation 𝑢′(𝑥⃗, 𝑡), such that:

𝑢(𝑥⃗, 𝑡) = 𝑢̅(𝑥⃗) + 𝑢′(𝑥⃗, 𝑡) (2.2)

The two components presented in Eq. 2.2 are shown in Fig. 2.2. It is easy
to show that:

𝑢̅ = 𝑢̅ 𝑢̅′ = 0 (2.3)

𝑢̅
𝑢′

Fig. 2.2.- Reynolds decomposition

Obviously, the same decomposition is valid for the other components of the
velocity and pressure:

:: 42 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑣 = 𝑣̅ + 𝑣′ 𝑤=𝑤
̅ + 𝑤′ 𝑝 = 𝑝̅ + 𝑝′ (2.4)

There is another kind of average, the ensemble average, which is more


appropriate for the analysis. In this case, the experiment is repeated many times
(𝑁), each repetition is called a realization. Thus, for realization 1, we have a
record 𝑢1 (𝑥⃗, 𝑡), for realization 2 we have a record 𝑢2 (𝑥⃗, 𝑡), and so on. An average
of the 𝑁 realizations for each (𝑥⃗, 𝑡) gives the ensemble average 〈𝑢(𝑥⃗, 𝑡)〉, defined
as:
𝑁
1
〈𝑢(𝑥⃗, 𝑡)〉 = ∑ 𝑢𝑖 (𝑥⃗, 𝑡) (2.2)
𝑁
𝑖=1
Note that the result of the ensemble average at a location 𝑥⃗, is not a value
as in the temporal average, but a function of time. This can be seen in Fig. 2.3.
For statistical processes that are ergodic and stationaries, ensemble and
temporal averages are equal:

〈𝑢〉 = 𝑢̅𝑖 (2.3)

In the analysis that follows we will consider that the experimental data
obtained in a turbulent flow behave statistically as an ergodic process strictly
stationary (i.e., not only the averages are the same but also any other statistical
property of the flow).
It is easy to show the following properties for any variable 𝑏 that is a
function of time
̅̅̅̅
𝜕𝑏 𝜕𝑏̅
= (2.4)
𝜕𝑥𝑖 𝜕𝑥𝑖

̅̅̅̅̅̅̅̅̅̅
∫ 𝑏 𝑑𝑥𝑖 = ∫ 𝑏̅ 𝑑𝑥𝑖 (2.5)

Replacing the Reynolds-decomposed variables (Eqs. 2.2 y 2.4) into


continuity equation (Eq. 87):

𝜕(𝑢̅ + 𝑢′) 𝜕(𝑣̅ + 𝑣′) 𝜕(𝑤


̅ + 𝑤′)
+ + =0 (2.6)
𝜕𝑥 𝜕𝑦 𝜕𝑧

:: 43 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

u
REALIZATION 1

𝑢1 (𝑡0 )

𝑡0 t
u REALIZATION 2

𝑢2 (𝑡0 )

𝑡0 t
. ⬚ .
. ⬚ .
. ⬚ .

u
REALIZATION 𝑁

𝑢𝑁 (𝑡0 )

𝑡0 t

ENSEMBLE
u AVERAGE
𝑁
1
〈𝑢(𝑡0 )〉 = ∑ 𝑢𝑖 (𝑡0 )
𝑁
𝑖=1

𝑡0 t

Fig. 2.3.- Ensemble average

Averaging Eq. 2.6:


̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕(𝑢̅ + 𝑢′) 𝜕(𝑣̅ + 𝑣′) 𝜕(𝑤 ̅ + 𝑤′)
+ + =0 (2.7)
𝜕𝑥 𝜕𝑦 𝜕𝑧

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕𝑢̅ 𝜕𝑣̅ 𝜕𝑤 ̅ 𝜕𝑢′ 𝜕𝑣′ 𝜕𝑤′
+ + + + + =0 (2.8)
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

:: 44 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕𝑢̅ 𝜕𝑣̅̅ 𝜕𝑤 ̅ 𝜕𝑤′


̅ 𝜕𝑢̅′ 𝜕𝑣′ ̅̅̅
+ + + + + =0 (2.9)
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

Using Eqs. 2.3:

𝜕𝑢̅ 𝜕𝑣̅ 𝜕𝑤
̅
+ + =0 (2.10)
𝜕𝑥 𝜕𝑦 𝜕𝑧

Substracting Eq. 2.10 to Eq. 2.6:

𝜕𝑢′ 𝜕𝑣′ 𝜕𝑤′


+ + =0 (2.11)
𝜕𝑥 𝜕𝑦 𝜕𝑧

Eqs. 2.10 and 2.11 indicate that the averaged velocities and their
fluctuations satisfy continuity.
Before repeating the same analysis with the Navier-Stokes equations, we
will modify slightly the equations. We will work first with the 𝑥 component of
the momentum equation. Let’s multiply the continuity equation by 𝑢:

𝜕𝑢 𝜕𝑣 𝜕𝑤
𝑢 +𝑢 +𝑢 =0 (2.12)
𝜕𝑥 𝜕𝑦 𝜕𝑧

Multiplying Eq. 2.12 by 𝜌 and adding it to the 𝑥 component of the Navier-


Stokes equation (Eq. 121):
𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑣 𝜕𝑤
𝜌( +𝑢 +𝑣 +𝑤 +𝑢 +𝑢 +𝑢 )
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧
𝜕𝑝 𝜕 2𝑢 𝜕 2𝑢 𝜕 2𝑢 (2.13)
=− + 𝜇 ( 2 + 2 + 2 ) + 𝜌𝑔𝑥
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧

We can recognize some product derivatives in the first parenthesis:

𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑢2
𝑢 +𝑢 = 2𝑢 =
𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑥
𝜕𝑢 𝜕𝑣 𝜕𝑢𝑣
𝑣 +𝑢 = (2.14)
𝜕𝑦 𝜕𝑦 𝜕𝑦
𝜕𝑢 𝜕𝑤 𝜕𝑢𝑤
𝑤 +𝑢 =
𝜕𝑧 𝜕𝑧 𝜕𝑧

:: 45 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Thus, the 𝑥 component of the Navier-Stokes is rewritten as


𝜕𝑢 𝜕𝑢2 𝜕𝑢𝑣 𝜕𝑢𝑤 𝜕𝑝 𝜕 2𝑢 𝜕 2𝑢 𝜕 2𝑢
𝜌( + + + )=− + 𝜇 ( 2 + 2 + 2 ) + 𝜌𝑔𝑥 (2.13)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧

Now, we will replace the Reynolds decomposed variables in Eq. 2.13: and
take the average. The resulting equation is:
̅̅̅̅̅̅̅̅̅̅̅
𝜕(𝑢̅ + 𝑢′) 𝜕(𝑢 ̅̅̅̅̅̅̅̅̅̅̅̅
̅ + 𝑢′)2 𝜕(𝑢 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅ + 𝑢′)(𝑣̅ + 𝑣′) 𝜕(𝑢 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅ + 𝑢′)(𝑤 ̅ + 𝑤′)
𝜌( + + + )
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧
̅̅̅̅̅̅̅̅̅̅
𝜕(𝑝̅ + 𝑝′) (2.14)
=− ̅̅̅̅̅̅̅̅̅̅̅
+ 𝜇∇2 (𝑢̅ + 𝑢′) + 𝜌𝑔𝑥
𝜕𝑥

The linear terms in Eq. 2.14 are easily decomposed into their average and
fluctuating parts:
̅̅̅̅̅̅̅̅̅̅̅
𝜕(𝑢 ̅ + 𝑢′) 𝜕𝑢̅̅ 𝜕𝑢′ ̅
= + =0 (2.15)
𝜕𝑡 𝜕𝑡 𝜕𝑡

̅̅̅̅̅̅̅̅̅̅
𝜕(𝑝̅ ̅ 𝜕𝑝̅
+ 𝑝′) 𝜕𝑝̅̅ 𝜕𝑝′
= + = (2.16)
𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑥

∇2 ̅̅̅̅̅̅̅̅̅̅̅ ̅ = ∇2 𝑢̅
(𝑢̅ + 𝑢′) = ∇2 𝑢̅ + ∇2 𝑢′ (2.17)

Let’s analyse now the non-linear terms:

̅̅̅̅̅̅̅̅̅̅̅̅
𝜕(𝑢̅ + 𝑢′)2 𝜕(𝑢 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
̅ 2 + 2𝑢̅𝑢′ + 𝑢′2 ) 𝜕𝑢 ̅̅̅ ̅̅̅̅̅̅
̅ 2 𝜕2𝑢 ̅̅̅̅
̅ 𝑢′ 𝜕𝑢′̅ 2 𝜕𝑢̅2 𝜕𝑢′
̅̅̅̅2
= = + + = + (2.18)
𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑥

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕(𝑢 ̅̅̅̅̅
̅ + 𝑢′)(𝑣̅ + 𝑣′) 𝜕𝑢̅𝑣̅ 𝜕𝑢′𝑣′
= + (2.19)
𝜕𝑦 𝜕𝑦 𝜕𝑦

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝜕(𝑢̅ + 𝑢′)(𝑤 ̅ + 𝑤′) 𝜕𝑢̅𝑤 ̅̅̅̅̅̅
̅ 𝜕𝑢′𝑤′
= + (2.20)
𝜕𝑧 𝜕𝑧 𝜕𝑧

Thus, the averaged 𝑥 component of the Reynolds equation becomes:

:: 46 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

̅̅̅̅2 𝜕𝑢̅𝑣̅ 𝜕𝑢′𝑣′


𝜕𝑢̅2 𝜕𝑢′ ̅̅̅̅̅ 𝜕𝑢̅𝑤 ̅̅̅̅̅̅
̅ 𝜕𝑢′𝑤′ 𝜕𝑝̅
𝜌( + + + + + )=− + 𝜇∇2 𝑢̅ + 𝜌𝑔𝑥 (2.21)
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑥

The right hand side of Eq. 2.21 can also be written as

𝜕𝑢̅2 𝜕𝑢̅𝑣̅ 𝜕𝑢̅𝑤


̅ 𝜕𝑢′ ̅̅̅̅2 𝜕𝑢′𝑣′
̅̅̅̅̅ 𝜕𝑢′𝑤′̅̅̅̅̅̅
+ + + + +
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧
𝜕𝑢̅ 𝜕𝑣̅ 𝜕𝑢̅ 𝜕𝑤̅ ̅̅̅̅2 𝜕𝑢′𝑣′
𝜕𝑢̅ 𝜕𝑢′ ̅̅̅̅̅
= 2𝑢̅ + 𝑢̅ + 𝑣̅ + 𝑢̅ +𝑤
̅ + + (2.22)
𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑥 𝜕𝑦
̅̅̅̅̅̅
𝜕𝑢′𝑤′
+
𝜕𝑧

𝜕𝑢̅2 𝜕𝑢̅𝑣̅ 𝜕𝑢̅𝑤


̅ 𝜕𝑢′ ̅̅̅̅2 𝜕𝑢′𝑣′
̅̅̅̅̅ 𝜕𝑢′𝑤′
̅̅̅̅̅̅
+ + + + +
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧
𝜕𝑢̅ 𝜕𝑣̅ 𝜕𝑤 ̅ 𝜕𝑢̅ 𝜕𝑢̅ ̅̅̅̅2 𝜕𝑢′𝑣′
𝜕𝑢̅ 𝜕𝑢′ ̅̅̅̅̅
= 𝑢̅ ( + + ) + 𝑢̅ + 𝑣̅ +𝑤
̅ + + (2.23)
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦
̅̅̅̅̅̅
𝜕𝑢′𝑤′
+
𝜕𝑧

Using Eq. 2.10, the term in parenthesis can be written as:

𝜕𝑢̅2 𝜕𝑢̅𝑣̅ 𝜕𝑢̅𝑤


̅ 𝜕𝑢′ ̅̅̅̅2 𝜕𝑢′𝑣′
̅̅̅̅̅ 𝜕𝑢′𝑤′
̅̅̅̅̅̅
+ + + + +
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧
̅̅̅̅2 𝜕𝑢′𝑣′
̅̅̅̅̅ 𝜕𝑢′𝑤′
̅̅̅̅̅̅ (2.24)
𝜕𝑢̅ 𝜕𝑢̅ 𝜕𝑢̅ 𝜕𝑢′
= 𝑢̅ + 𝑣̅ +𝑤 ̅ + + +
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

Eq. 2.21 becomes:

𝜕𝑢̅ 𝜕𝑢̅ ̅̅̅̅2 𝜕𝑢′𝑣′


𝜕𝑢̅ 𝜕𝑢′ ̅̅̅̅̅ 𝜕𝑢′𝑤′
̅̅̅̅̅̅ 𝜕𝑝̅
𝜌 (𝑢̅ + 𝑣̅ +𝑤
̅ + + + )=− + 𝜇∇2 𝑢̅ + 𝜌𝑔𝑥 (2.25)
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥

Remember that Eq. 2.21 is not more that Newton’s second law applied to
an incompressible Newtonian fluid and expressed in terms of temporal mean
values. In order to interpret the meaning of Eq. 2.21, we should remember the
equation of Newton’s second law: 𝑑(𝑚𝑣⃗)⁄𝑑𝑡 = 𝐹⃗ . As the mass is preserved, we
can write 𝑚𝑑(𝑣⃗)⁄𝑑𝑡 = 𝐹⃗ . For simplicity, let’s consider the 𝑥 component:

:: 47 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑑𝑢
𝑚 = 𝐹𝑥 (2.26)
𝑑𝑡

Of course, Newton’s second law is not applied to fluids in the form of Eq.
2.26, but as Eq. 51 (integral approach) or Eq. 124 (differential approach), but it
is used in this explanation for the sake of clarity. The problem with Eq. 2.26 (and
the reason why in fluids is expressed in a different way is to define the mass 𝑚
in a flow. Anyway, we can divide Eq. 2.26 by a volume and work with forces per
unit volume, 𝐹𝑉𝑥 ≡ 𝐹𝑥 ⁄𝑉 :

𝑑𝑢
𝜌 = 𝐹𝑉𝑥 (2.27)
𝑑𝑡

We can identify 𝐹𝑉𝑥 in Eq. 124 as the resulting of the force due to pressure,
viscosity, and gravity. We would like to work with a temporal mean acceleration,
and write Eq. 2.27 as:

𝑑𝑢̅
𝜌 = 𝐹̅𝑥 (2.28)
𝑑𝑡

We can define an acceleration based on the temporal mean quantities. For


a statistically stationary flow this acceleration is:

𝑑𝑢̅ 𝜕𝑢̅ 𝜕𝑢̅ 𝜕𝑢̅


= 𝑢̅ + 𝑣̅ +𝑤
̅ (2.29)
𝑑𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧

Thus, introducing Eq. 2.29 into Eq. 2.25:


̅̅̅̅2 𝜕𝑢′𝑣′
𝑑𝑢̅ 𝜕𝑢′ ̅̅̅̅̅ 𝜕𝑢′𝑤′
̅̅̅̅̅̅ 𝜕𝑝̅
𝜌( + + + )=− + 𝜇∇2 𝑢̅ + 𝜌𝑔𝑥 (2.30)
𝑑𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥

We can see that Eq. 2.30 is very similar to Eq. 2.28. The only problem is
the terms due to the fluctuations in the parenthesis. In order to have only the
variation of momentum in the left hand side of Eq. 2.30, we take the terms due
to the fluctuations towards the right hand side of the equation:

𝑑𝑢̅ 𝜕𝑝̅ ̅̅̅̅2 𝜕𝑢′𝑣′


𝜕𝑢′ ̅̅̅̅̅ 𝜕𝑢′𝑤′
̅̅̅̅̅̅
𝜌 =− + 𝜇∇2 𝑢̅ + 𝜌𝑔𝑥 − 𝜌 ( + + ) (2.31)
𝑑𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧

To pass the terms due to the fluctuations from one side of the equal
sign to the other side is much more than an algebraic step. It changes the
interpretation that we can give to the average of products of the fluctuations. As

:: 48 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

they are now in the right side of the equation, we can interpret them as forces
arising from the turbulence. Thus, the terms −𝜌𝑢 ̅̅̅̅
′2 , −𝜌𝑢 ′ 𝑣 ′ and −𝜌𝑢
̅̅̅̅̅̅ ′ 𝑤 ′ are
̅̅̅̅̅̅
called Reynold’s apparent stresses or, simply, Reynold’s stresses , or turbulent
stresses. Let’s work with the terms associated to the viscous stresses 𝜇∇2 𝑢̅ :

𝜕 2 𝑢̅ 𝜕 2 𝑢̅ 𝜕 2 𝑢̅
∇2 𝑢̅ = + + (2.32)
𝜕𝑥 2 𝜕𝑦 2 𝜕𝑧 2

As the continuity equation is ∇ ∙ 𝑣⃗ = 0,, we can take its derivative


with respect to 𝑥 and add to Eq. 2.32, and it will not change, i.e.:

𝜕
∇2 𝑢̅ = ∇2 𝑢̅ + ∇ ∙ 𝑣̅⃗ (2.33)
𝜕𝑥

𝜕 2 𝑢̅ 𝜕 2 𝑢̅ 𝜕 2 𝑢̅ 𝜕 2 𝑢̅ 𝜕 2 𝑣̅ 𝜕 2𝑤
̅
∇2 𝑢̅ = 2
+ 2
+ 2
+ + + (2.34)
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥𝜕𝑥 𝜕𝑥𝜕𝑦 𝜕𝑥𝜕𝑧

𝜕 𝜕𝑢̅ 𝜕 𝜕𝑢̅ 𝜕 𝜕𝑢̅ 𝜕 𝜕𝑢̅ 𝜕 𝜕𝑣̅ 𝜕 𝜕𝑤


̅
∇2 𝑢̅ = + + + + + (2.35)
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑧 𝜕𝑥

𝜕 𝜕𝑢̅ 𝜕 𝜕𝑢̅ 𝜕𝑣̅ 𝜕 𝜕𝑢̅ 𝜕𝑤


̅
∇2 𝑢̅ = (2 ) + ( + )+ ( + ) (2.36)
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑧 𝜕𝑧 𝜕𝑥

Multiplying Eq. 2.36 by 𝜇:

𝜕 𝜕𝑢̅ 𝜕 𝜕𝑢̅ 𝜕𝑣̅ 𝜕 𝜕𝑢̅ 𝜕𝑤


̅
𝜇∇2 𝑢̅ = (2𝜇 ) + (𝜇 ( + )) + (𝜇 ( + )) (2.37)
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑧 𝜕𝑧 𝜕𝑥

We identify the terms in parenthesis as some of the viscous stress


that appear in Eq. 114. Denoting with the sub-index 𝑉 to denote “viscous”, we
have:

𝜕𝑢̅ 𝜕𝑢̅ 𝜕𝑣̅ 𝜕𝑢̅ 𝜕𝑤


̅
𝜏𝑉𝑥𝑥 = 2𝜇 , 𝜏𝑉𝑦𝑥 = 𝜇 ( + ) , 𝜏𝑉𝑧𝑥 = 𝜇 ( + ) (2.38)
𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑧 𝜕𝑥

:: 49 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

We can write now:

𝜕𝜏𝑉𝑥𝑥 𝜕𝜏𝑉𝑦𝑥 𝜕𝜏𝑉𝑧𝑥


𝜇∇2 𝑢̅ = + + (2.39)
𝜕𝑥 𝜕𝑦 𝜕𝑧

Replacing in Eq. 2.31:

𝑑𝑢̅ 𝜕𝑝̅ 𝜕𝜏𝑉𝑥𝑥 𝜕𝜏𝑉𝑦𝑥 𝜕𝜏𝑉𝑧𝑥 ̅̅̅̅2 𝜕𝑢′𝑣′


𝜕𝑢′ ̅̅̅̅̅ 𝜕𝑢′𝑤′
̅̅̅̅̅̅
𝜌 =− + + + + 𝜌𝑔𝑥 − 𝜌 ( + + ) (2.40)
𝑑𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝑑𝑢̅ 𝜕𝑝̅ 𝜕 ̅̅̅̅2 ) + 𝜕 (𝜏𝑉𝑦𝑥 − 𝜌𝑢′𝑣′


̅̅̅̅̅) + 𝜕 (𝜏𝑉𝑧𝑥 − 𝜌𝑢′𝑤′
̅̅̅̅̅̅)
𝜌 =− + (𝜏𝑉𝑥𝑥 − 𝜌𝑢′
𝑑𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧 (2.41)
+ 𝜌𝑔𝑥

Using the sub-index 𝑇 to denote the turbulent or Reynolds stresses:

̅̅̅̅2 ,
𝜏 𝑇𝑥𝑥 = −𝜌𝑢′ ̅̅̅̅̅ ,
𝜏 𝑇𝑦𝑥 = −𝜌𝑢′𝑣′ ̅̅̅̅̅̅
𝜏 𝑇𝑧𝑥 = −𝜌𝑢′𝑤′ (2.42)

We define the total stress as the sum of the viscous and the turbulent one:

𝑇𝑥𝑥 = 𝜏𝑉𝑥𝑥 + 𝜏 𝑇𝑥𝑥 , 𝑇𝑦𝑥 = 𝜏𝑉𝑦𝑥 + 𝜏 𝑇𝑦𝑥 , 𝑇𝑧𝑥 = 𝜏𝑉𝑧𝑥 + 𝜏 𝑇𝑧𝑥 (2.43)

In general:

𝑇𝑖𝑗 = 𝜏𝑉𝑖𝑗 + 𝜏 𝑇𝑖𝑗 (2.44)

Where the viscous stress (associated to the temporal mean velocities) is


given by:

𝜕𝑢̅𝑖 𝜕𝑢̅𝑗
𝜏𝑉𝑖𝑗 = 𝜇 ( + ) (2.45)
𝜕𝑥𝑗 𝜕𝑥𝑖

And the turbulent stress is:

̅̅̅̅̅̅
𝜏 𝑇𝑖𝑗 = −𝜌𝑢 ′ ′
𝑖 𝑢𝑗 (2.46)

Thus, the momentum equation in the 𝑥direction for a stationary turbulent


flow, in terms of the temporal mean values is given by:

:: 50 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕𝑢̅ 𝜕𝑢̅ 𝜕𝑢̅ 𝜕𝑝̅ 𝜕𝑇𝑥𝑥 𝜕𝑇𝑦𝑥 𝜕𝑇𝑧𝑥


𝜌 (𝑢̅ + 𝑣̅ +𝑤
̅ )=− + + + + 𝜌𝑔𝑥 (2.47)
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧

In the same way, we can obtain the 𝑦 and 𝑧 components:


𝜕𝑣̅ 𝜕𝑣̅ 𝜕𝑣̅ 𝜕𝑝̅ 𝜕𝑇𝑥𝑦 𝜕𝑇𝑦𝑦 𝜕𝑇𝑧𝑦
𝜌 (𝑢̅ + 𝑣̅ +𝑤
̅ )=− + + + + 𝜌𝑔𝑦 (2.48)
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕𝑤
̅ 𝜕𝑤
̅ 𝜕𝑤
̅ 𝜕𝑝̅ 𝜕𝑇𝑥𝑧 𝜕𝑇𝑦𝑧 𝜕𝑇𝑧𝑧
𝜌 (𝑢̅ + 𝑣̅ +𝑤
̅ )=− + + + + 𝜌𝑔𝑧 (2.49)
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

Eqs. 2.47, 2.48 and 2.49 are the Reynolds equations for turbulent flows,
and they constitute an important advance in the study and analysis of the
turbulence. In his paper, Reynolds also derived and discussed the equation for
the turbulent kinetic energy.
Although we have had some progress in the analysis of turbulent flows,
we are not in conditions to solve any problem yet. We have a system formed by
4 partial differential equations (with their boundary conditions): Eqs. 2.10, 2.47,
̅, 𝑝̅ , ̅̅̅̅
2.48 and 2.49. But we have 10 unknowns: 𝑢̅, 𝑣̅ , 𝑤 𝑢′2 , ̅̅̅̅
𝑣′2 , ̅̅̅̅̅
𝑤′2 , ̅̅̅̅̅
𝑢′𝑣′ , ̅̅̅̅̅̅
𝑢′𝑤′ , ̅̅̅̅̅̅
𝑣′𝑤′. We
cannot solve any problem of turbulence using the equations derived with the
Reynolds decomposition approach if we do not know relations for −𝜌𝑢 𝑖 𝑢𝑗 . This is
̅̅̅̅̅̅
′ ′

the so-called “closure problem of turbulence”. Basically, the problem is how to


model the turbulent stresses. As there is not a theory based only on the first
principles of the physics, all the available models necessarily require some
experimental data. We will present two closures of the problem, both of them
widely used in engineering. They are the Boussinesq’s eddy viscosity model and
Prandtl mixing length theory.
Boussinesq closure of the turbulence: Eddy viscosity
Boussinesq published in 1877 a compilation of his research on water
flows. (The book can be read in the site
http://gallica.bnf.fr/ark:/12148/bpt6k56673076 , but cannot be downloaded). In
Eqs. 12 of his book, Boussinesq presents the stresses in a similar form (but with
slightly different notation) than those given in Eq. 114 in tensorial notation, or
Eqs. 115 to 120 (see Fig. 2.4). The big difference is that in the later equations,
appears 𝜇, and in Boussinesq’s equations there is a coefficient 𝜀. But before to
discuss the meaning and value of 𝜀, it is interesting to note that Boussinesq was
solving a problem that arose almost 20 years later, with Reynolds. It is
interesting to note that Boussinesq already considered temporal averages, but

:: 51 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

his mistake was to assume that the velocity fluctuations were not correlated, i.e.,
𝑢𝑖′ 𝑢𝑗′ = 0 (using our notation). In this way, he lost additional components to the
̅̅̅̅̅̅
stresses. However, it was known that the equations of Navier-Stokes provided
results in agreement with the experiments for conduits with small flow area, but
it failed for larger conduits. In the latter case, additional effects appear in the
flow, with the final effect that the viscosity seems to be larger (in our language
we would say that the Navier-Stokes were in agreement with measurements in
laminar flows, but other effects should be taking into account when dealing with
turbulent flows. Details are in the first 50 pages of Boussinesq’s book. For open
channel flows he gives:

𝜀 = 𝜌𝑔𝐴ℎ𝑢0 (2.50)

Eq. 2.50 corresponds to Boussinesq’s Eq. 13. In an open channel flow, ℎ is the
flow depth and 𝑢0 is the “velocity at the wall”. 𝐴 is a coefficient that depends on
the wall roughness varies little with ℎ and 𝑢0 . Note that in order that Eq. 2.50
be dimensionally homogeneous, the dimensions of 𝐴 must be T2L-1
(𝑡𝑖𝑚𝑒 𝑙𝑒𝑛𝑔𝑡ℎ).
2⁄

Fig. 2.4.- Equations for the stresses presented by Boussinesq


in 1877 in his Essai sur la théorie des eaux courantes. 𝜺 is the
eddy viscosity.

In our current language, 𝜀 is named eddy viscosity, or turbulent viscosity.


Thus, according to Boussinesq, the turbulent stresses 𝜏 𝑇𝑖𝑗 , can be computed in
the same way than the viscous stresses, using the eddy viscosity 𝜀 instead of 𝜇,
the molecular dynamic viscosity. The turbulent stresses become:

𝜕𝑢̅
𝜏 𝑇𝑥𝑥 = 𝜀2 (2.51)
𝜕𝑥

𝜕𝑣̅
𝜏 𝑇𝑦𝑦 = 𝜀2 (2.52)
𝜕𝑦

:: 52 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕𝑤
̅
𝜏 𝑇𝑧𝑧 = 𝜀2 (2.53)
𝜕𝑧

𝜕𝑢̅ 𝜕𝑣̅
𝜏 𝑇𝑥𝑦 = 𝜏 𝑇𝑦𝑥 = 𝜀 ( + ) (2.54)
𝜕𝑦 𝜕𝑥

𝜕𝑢̅ 𝜕𝑤
̅
𝜏 𝑇𝑥𝑧 = 𝜏 𝑇𝑧𝑥 = 𝜀 ( + ) (2.55)
𝜕𝑧 𝜕𝑥

𝜕𝑣̅ 𝜕𝑤
̅
𝜏 𝑇𝑦𝑧 = 𝜏 𝑇𝑧𝑦 = 𝜀 ( + ) (2.60)
𝜕𝑧 𝜕𝑦

There is a strong difference between 𝜇 and 𝜀: the dynamic viscosity 𝜇 is a property


of the fluid and the eddy viscosity 𝜀 is a property of the flow. Thus, if we are
working with water at a given temperature, we can look for the value of the
viscosity in any book and use it, independently if the flow is in a cylindrical or
square pipe, a rectangular or trapezoidal channel. On the contrary, the eddy
viscosity depends on the flow and the geometry and it is independent of the fluid
(at least, for flows with high Reynolds numbers). As it is was said in the
Introduction, one of the features of the turbulent flows is they are highly efficient
in the momentum transfer process, when compared with laminar flows. Thus,
another characteristic of the eddy viscosity is that it is much larger (several
orders of magnitude) than the dynamic viscosity:𝜀 ≫ 𝜇.
In analogy to the molecular kinematic viscosity  = 𝜇 ⁄𝜌, a kinematic eddy
viscosity is also used  𝑇 = 𝜀 ⁄𝜌 . Eddy viscosities for some flow configurations are
presented in Table 2.1. Known the eddy viscosity, we can solve the problem of
getting the velocity distribution in a turbulent flow. We have to solve the system
of partial differential equations formed by continuity and Reynolds equations, in
addition to an expression for the eddy viscosity.
WARNING! In the current technical literature 𝜀 is frequently used to design the
kinematic eddy viscosity (i.e.,  𝑇 = 𝜀) and not the dynamic turbulent viscosity, as
Boussinesq designed in his original work and shown in Fig. 2.4..
As an example of Boussinesq’s eddy viscosity approach to close Eqs. 2.47
to 2.49, we can solve the problem of the 2-D steady, uniform, turbulent flow on
an inclined plane that we have been using as example along this notes (Fig. 18).
Using the same arguments than in the laminar case, the equations of continuity
(Eq. 2.10) and momentum (Reynolds equations, Eqs. 2.47 to 2.49) are reduced to:

𝜕𝑢̅
=0 (2.61)
𝜕𝑥

:: 53 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Reynolds equation in the 𝑥 direction:


𝜕𝑇𝑦𝑥
0= + 𝜌𝑔𝑥 (2.62)
𝜕𝑦

Reynolds equation in the 𝑦 direction:

Table 2.1.- EDDY VISCOSITY FOR SOME FLOW CONFIGURATIONS

FLOW CONFIGURATION EDDY VISCOSITY

2-D Open channel flow

 𝑇 = 𝐻𝑢∗ (1 − )
𝑦
 = 0.4 ;  =
𝐻

Axisymetric jet

 𝑇 = 0.013𝑉0 𝑑0

Central region of the flow in a pipe


𝑇 𝐶 𝑓
= 𝑅𝑒√
 2 8
𝑈𝐷
𝐶 ≈ 0.07 ; 𝑅𝑒 =

𝜕𝑝̅
0=− + 𝜌𝑔𝑦 (2.63)
𝜕𝑦

From where we obtain that the mean pressure distribution is hydrostatic.


With the boundary condition 𝑦 = 𝐻 , 𝑝̅ = 0 (relative pressure), we obtain:

:: 54 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑝̅ = 𝜌𝑔 cos 𝜃 (𝐻 − 𝑦) (2.64)

Reynolds equation along 𝑥 (Eq. 2.62) can be integrated once with respect
to 𝑦:

𝑇𝑦𝑥 = −𝜌𝑔𝑥 𝑦 + 𝐶1 (2.65)

As there is not shear acting on the free surface, the boundary condition at
𝑦 = 𝐻 is 𝑇𝑦𝑥 = 0. Thus:

𝐶1 = 𝜌𝑔𝑥 𝐻 (2.66)

Using Eq. 2.43 and the value of 𝐶1 in Eq. 2.65:

𝜏𝑉𝑦𝑥 + 𝜏 𝑇𝑦𝑥 = 𝜌𝑔𝑥 (𝐻 − 𝑦) (2.67)

Replacing the shear stresses by Eqs. 118 and 2.54:


𝜕𝑢̅
(𝜇 + 𝜀) = 𝜌𝑔𝑥 (𝐻 − 𝑦) (2.68)
𝜕𝑦

Considering that the molecular viscosity 𝜇 is much smaller than the eddy
viscosity 𝜀:
𝜕𝑢̅
𝜀 = 𝜌𝑔𝑥 (𝐻 − 𝑦) (2.69)
𝜕𝑦

To drop 𝜇 from the momentum equation limits its use only to the region
where turbulence dominates over viscosity. We will see later that even in highly
turbulent flows bounded by smooth walls, in a small region near the wall
molecular viscosity effects are important.
In order to integrate Eq. 2.69, we have to decide what value of 𝜀 will be
used. We can work with the expression given by Boussinesq (Eq. 2.50) or that
presented in Table 2.1.
Solution considering Boussinesq expression for 𝜀:
Eq. 2.50 is 𝜀 = 𝜌𝑔𝐴ℎ𝑢0 , and Eq. 2.69 becomes:
𝜕𝑢̅ 𝑔𝑥
= (𝐻 − 𝑦) (2.70)
𝜕𝑦 𝜌𝑔𝐴ℎ𝑢0

:: 55 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Integrating:
𝑔𝑥 1
𝑢̅ = (𝐻𝑦 − 𝑦 2 ) + 𝐶2 (2.71)
𝑔𝐴ℎ𝑢0 2

Although still have to determine the integration constant 𝐶2 , the


velocity distribution for a turbulent flow given by Eq. 2.71 is essentially the
same than the distribution obtained for a laminar flow (Eq. 149) (parabolic
distribution). This should not surprise us: we only changed the value of the
constant that multiplies the derivative of the velocity. Determination of 𝐶2 is
conceptually more complicated. We cannot impose the condition that at 𝑦 = 0,
𝑢̅ = 0 because very close to the bottom Eq. 2.69 is not valid. We should know the
distance 𝑦𝑇 from the wall where turbulence dominates, and the value of the
velocity, 𝑢𝑇 , at that location. Instead of doing that, we will follow Boussinesq
ideas. According to him, in there is a velocity at the wall, 𝑢0 , which is rather
large. In a footnote of his book (p. 51) he works a typical value of 1 m/s. Thus,
according to Boussinesq, at 𝑦 = 0, 𝑢̅ = 𝑢0 , and the resulting velocity distribution
is:
sin 𝜃 1
𝑢̅ = 𝑢0 + (𝐻𝑦 − 𝑦 2 ) (2.72)
𝐴ℎ𝑢0 2

A comparison of the laminar and Boussinesq’s turbulent velocity profiles are


sketched in Fig. 2.4. A feature of the turbulent velocity profiles is that they are
flatter than the laminar ones. In Boussinesq constant eddy viscosity model this
is achieved by the term 𝑢0 , the fluid slip at the wall.

Fig. 2.4.- Comparison of the velocity distribution for laminar flow and
turbulent flow according to the eddy viscosity model of Boussinesq
(1877) that includes a velocity 𝒖𝟎 at the wall.

:: 56 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Solution considering the expression given in Table 2.1for 𝜀:

In this case, the eddy viscosity is given by:


𝑦
𝜀 = 𝜌𝑢∗ (𝐻 − 𝑦) (2.73)
𝐻

The equation of motion gives:


𝜕𝑢̅ 𝜌𝑔𝑥 𝑔𝑥 𝐻
= (𝐻 − 𝑦) = (2.74)
𝜕𝑦 𝜀 𝑢∗ 𝑦

Integrating:
𝑔𝑥 𝐻
𝑢̅ = ln(𝑦) + 𝐶3 (2.75)
𝑢∗

Again, we cannot use a boundary condition at 𝑦 = 0 to determine the


integration constant 𝐶3 . Mathematically, the logarithm in Eq. 2.75 blows up.
Physically, the velocity distribution is valid only in the flow region where
turbulence dominates (we have neglected 𝜇). 𝐶3 is determined fitting Eq. 2.75 to
experimental data. Using Eqs. 140 and 156 together to the definition of shear
velocity (Eq. 67) 𝑔𝑥 𝐻 = 𝑔𝐻 sin 𝜃 = 𝜏0 ⁄𝜌 = 𝑢∗2 :
𝑢∗
𝑢̅ = ln(𝑦) + 𝐶3 (2.76)

The logarithmic velocity distribution given by Eq. 2.76 is sketched in
Fig.2.5. It is observed that, for 𝑦 smaller than a certain value, the velocity
becomes negative and 𝑢̅ → −∞ when 𝑦 → 0. Obviously, the velocity distribution
is valid only in the region where the flow is turbulent, i.e., for 𝑦 farther than a
specific distance from the bottom.
It is worth to insist that, the velocity profiles obtained for turbulent flow
relay on experimental data. The parabolic profile of Boussinesq requires to know
experimentally the coefficient 𝐴 and to have a method to determine 𝑢0 . The
logarithmic profile require experiments to determine  and 𝐶3 . The opposite
happens with the velocity profile for laminar flow. The parabolic profile obtained
from the Navier-Stokes equation and given by Eq. 149 only requires to know the
properties of the fluid, the flow and the acceleration due to gravity.

:: 57 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Fig. 2.5.- Logarithmic velocity profile. The velocity distribution is


valid from a certain distance of the bottom, where the flow is
turbulent.

Prandtl closure of the turbulence: Mixing length


Ludwig Prandtl is one of the great researchers in fluid mechanics of the
XX century. His mixing length concept (or theory), is just one of his many
contributions to the study of turbulence. It was proposed in 1925 (Prandtl, 1925).
In the mixing length theory, Prandtl assumes that portions of fluid can
move due to the velocity fluctuations a distance 𝑙 without losing its identity,
basically preserving its momentum (per unit volume). To fix ideas, let’s consider
a 2-D flow with a mean velocity distribution as that sketched in Fig.2.6. In that

𝑢̅(𝑦) 𝑢′(𝑦)
𝑦+𝑙
𝑢̅(𝑦 + 𝑙) B

𝑣′(𝑦)

𝑦
𝑢̅(𝑦) A

𝑢̅
Fig. 2.6.- Concept of mixing length
figure, the continuous curve represents the mean velocity profile and the circles

:: 58 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

a mass of fluid. A mass of fluid that initially was at A with velocicty 𝑢̅(𝑦), due to
the vertical fluctuation 𝑣′ is transported to the location B, without losing its
identity, i.e., preserving its momentum in the 𝑥 direction (per unit volume) 𝜌𝑢̅.
The distance of the displacement of the mass of fluid is 𝑙, named mixing length
by Prandtl. As the mass of fluid did not changed its properties during the
displacement along the 𝑦 axis, when it arrives to the location B, has the velocity
𝑢̅(𝑦), which is imposed at the level (𝑦 + 𝑙). It is easy to see that the difference
between this new velocity at (𝑦 + 𝑙 ) and the mean velocity at this level, 𝑢̅(𝑦 + 𝑙),
corresponds to the velocity fluctuation 𝑢′:

𝑢′ = 𝑢̅(𝑦) − 𝑢̅(𝑦 + 𝑙) (2.77)

Expanding in Taylor’s series and neglecting terms of second order and


superior:
𝜕𝑢̅
𝑢′ = 𝑢̅(𝑦) − (𝑢̅(𝑦) + 𝑙 +⋯) (2.78)
𝜕𝑦

Thus:
𝜕𝑢̅
𝑢′ = −𝑙 (2.79)
𝜕𝑦

Note that, for the mean velocity profile of Fig. 2.6, 𝜕𝑢̅⁄𝜕𝑦 > 0, and 𝑢′ < 0 .
The fluctuating velocity in the 𝑦 direction was positive, 𝑣 ′ > 0. This means that

𝑢′ 𝑣 ′ < 0 (2.80)

Experimental measurements show that 𝑢′ and 𝑣 ′ are well correlated, i.e.,


they are of the same order of magnitude:

|𝑢′ |~|𝑣 ′ | (2.81)

From the las two arguments (Eqs. 2.80 and 2.81), we can write:
𝜕𝑢̅
𝑣 ′ ~𝑙 (2.82)
𝜕𝑦

̅̅̅̅̅
Thus, the turbulent shear stress 𝜏 𝑇𝑥𝑦 = −𝜌𝑢 ′ 𝑣′ (Eq. 2.42) can be written

as:

:: 59 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝜕𝑢̅ 2
𝜏 𝑇𝑥𝑦 2
= 𝜌𝑙 ( ) (2.83)
𝜕𝑦

From Eq. 2.83 we obtain that the eddy viscosity for this flow is given by:
𝜕𝑢̅
𝜀 = 𝜌𝑙 2 (2.84)
𝜕𝑦

We can try to solve the problem of the permanent uniform flow over an
inclined plane using Prandtl’s expression for the turbulent shear stress (Eq.
2.83). We had obtained (Eq. 2.67):

𝜏𝑉𝑦𝑥 + 𝜏 𝑇𝑦𝑥 = 𝜌𝑔𝑥 (𝐻 − 𝑦) (2.67)

Eqs. 2.38 and 2.83:

𝜕𝑢̅ 𝜕𝑢̅ 2
𝜇 2
+ 𝜌𝑙 ( ) = 𝜌𝑔𝑥 (𝐻 − 𝑦) (2.85)
𝜕𝑦 𝜕𝑦

If we consider the region of the flow where turbulence dominates:

𝜕𝑢̅ 2
𝑙2 ( ) = 𝑔𝑥 (𝐻 − 𝑦) (2.86)
𝜕𝑦

We cannot go further than Eq. 2.86 if we do not know an expression for


the mixing length 𝑙. At this point we need to know 𝑙. Although the story is not
simple, and involves the competition between Pradtl’s research group in
Gottingen and one of his former Ph.D. students, emigrated to the USA, Theodor
von Kármán. In a paper presented in 1930 he assumes self-similarity of the
velocity profiles and proposes (Karman, 1930):
𝑑𝑢̅
𝑑𝑦
𝑙=  2 (2.87)
𝑑 𝑢̅
𝑑𝑦 2
The proportionality coefficient  is named von Karman coefficient and
must be determined experimentally. Von Karman recognized that Eq. 2.87 is not
valid near the walls, because in those regions the turbulence is damped and
viscosity needs to be taking into account. Von Karman applies his theory to the
flow between two parallel plates separated a distance 2𝐻, as the centre-line is a
symmetry axis, the problem is equivalent to the flow over an inclined plane.
Replacing Eq. 2.87 in Eq. 2.86:

:: 60 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑑𝑢̅ 4
( )
𝑑𝑦
2 2 = 𝑔𝑥 (𝐻 − 𝑦) (2.88)
𝑑 2 𝑢̅
( 2)
𝑑𝑦

Making the change of variables 𝑦 ′ = 𝐻 − 𝑦 and calling 𝑈′ = 𝑑𝑢̅⁄𝑑𝑦 =


− 𝑑𝑢̅⁄𝑑𝑦′ it is easy to perform the integration. Without going into details, and
copying the result from von Karman’s paper:
𝑑𝑢̅ √𝑔𝑥 1
𝑈′ = = (2.89)
𝑑𝑦′ 2 √𝐻 − √𝑦′

A word regarding the integration constant needs to be mentioned at this


point. Following to von Karman, when we approach to the wall 𝑑𝑢̅⁄𝑑𝑦 becomes
very large because 𝜇 is very small. Thus, he computes the integration constant
evaluating at 𝑦 = 0 , (𝑦 ′ = 𝐻) and impossing 𝑈′ → ∞. Eq. 2.89 is integrated again,
with the boundary condition that the velocity is maximum (𝑈𝑚𝑎𝑥 ) at t 𝑦 = 𝐻 ,
(𝑦 ′ = 𝐻). Using 𝑢∗ = √𝑔𝐻 sin 𝜃 the velocity distribution is:

𝑢∗ 𝑦 𝑦
𝑢̅ = 𝑈𝑚𝑎𝑥 + [ln (1 − √1 − ) + √1 − ] (2.90)
 𝐻 𝐻

Note that Eq. 2.90 fails at 𝑦 = 0. This is correct because it is valid only in
the region where the turbulent stresses dominate over the viscous ones.
Further, von Karman assumes that 𝑙 should “quietly diminish to zero and
the velocity distribution near to the wall becomes” (Eq. 25 in von Karman’s
paper):

𝜕𝑢̅ 𝜕𝑢̅ 2
𝜏𝑦𝑥 − 𝜇 = 𝜌(𝑦)2 ( ) (2.91)
𝜕𝑦 𝜕𝑦

Von Karman recognizes that a Japanese researcher, Wada, got the


equation in 1927, but the article was practically unknown because it was
published in Japanese in a journal of naval architecture. According to von
Karman, Wada applied Eq. 2.91 to the whole channel “which made the formulas
a little too complicated”. Von Karman’s results are already contained in Wada’s
work “in an implied form”. The right hand side term of Eq. 2.91 is the turbulent
shear stress. A comparison with Eq. 2.83, gives the mixing length in the
turbulent region near the wall (or bottom):

:: 61 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑙 = 𝑦 (2.92)

In order to clarify the applicability of Eq. 2.92, consider the flow on a


smooth wall, as sketched in Fig. 2.7. The existence of a “laminar layer” in contact
with a smooth wall was known since the 1920’s (von Karman, 1930). Due to the
presence of the wall, the turbulent fluctuations are damped close to the wall,
generating a region in which the viscous effects dominate over the turbulent ones
(𝜏𝑉𝑖𝑗 ≪ 𝜏 𝑇𝑖𝑗 ). Now we call that region viscous sublayer. Far enough of the wall,
turbulence dominates over viscosity and we have the turbulent region, where
𝜏 𝑇𝑖𝑗 ≪ 𝜏𝑉𝑖𝑗 . The passage from the viscous sublayer to the turbulent region is not
abrupt and a buffer layer is identified in which viscous and turbulent effects are
equally important (𝜏𝑉𝑖𝑗 ~𝜏 𝑇𝑖𝑗 ). Further, the turbulent region can be divided in an
inner region where there is a local energy equilibrium, and an outer region, which
is dominated by large eddies which transport low momentum and high turbulent
energy from the outer part of the inner region towards the upper part of the flow
(Bradshaw, 1972).

OUTER TURBULENT
REGION REGION
FLOW

𝑦 INNER
REGION
BUFFER LAYER
VISCOUS SUBLAYER

Fig. 2.7.- Regions of the turbulent flow over a smooth bottom.

The applicability of Eq. 2.92 is near the wall, but in the turbulent region.
Using 𝑙 = 𝑦 in Eq. 2.86, we have:

𝜕𝑢̅ 2
(𝑦) ( ) = 𝑔𝑥 (𝐻 − 𝑦)
2 (2.93)
𝜕𝑦
At this point if the analysis, we should mention that according to
Schlichting (1979, p. 587) Eq. 2.92 was Prandtl’s idea that was presented in a
couple of papers in 1925 and 1926 (As I do not read German, I cannot verify
Schlichting’s statement, but both papers do not contain an explicit equation
similar to Eq. 2.92. It is possible that a linearly decreasing mixing length near
the wall be mentioned in the text. Neither in von Karman's article there is an

:: 62 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

explicit expression for the mixing length, but it is inferred from his Eq. 25, as
indicated before).
Limiting the application of Eq. 2.93 to a region very near the wall, 𝑦 ≪ 𝐻,
and recalling that 𝜏0 = 𝜌𝑔𝐻 sin 𝜃, Eq. 2.93 is simplified to:
𝑑𝑢̅ 𝜏0
 = √ (2.94)
𝑑𝑦 𝜌

According to Schlichting (1979, p.588) to considerer that the shear stress


was constant and equal to that existing at the bottom was a “far-reaching
assumption” introduced by Prandtl.
Integrating Eq. 2.94 and using the definition of shear velocity, the velocity
distribution is:
𝑢∗
𝑢̅ = ln(𝑦) + 𝐶 (2.95)

Eq. 2.95 is named the logarithmic law for the velocity distribution or,
simply, logarithmic law. The constants and 𝐶 needs to be determined from
experimental data. We should not be surprised to see that we got the same result
as that obtained with the eddy viscosity. The value of 𝜀 used to obtain the
logarithmic profile given by Eq. 2.76 was obtained using the solution resulting
from the mixing length (Eq. 2.95).
Although there is a discussion about the universality of von Karman’s
constant, it is usually taken as  = 0.4. A dependency on the suspended
sediments concentration exists. Determination of the constant 𝐶 deserves a
separated analysis, because it depends on the characteristics of the wall.
In addition to the mixing length given by Eq. 2.92, other have been
proposed. For example, for the turbulent flow in a pipe of radius 𝑅, Nikuradse
suggested (Schlichting, 1979):

𝑙 𝑦 2 𝑦 4
= 0.14 − 0,08 (1 − ) − 0,06 (1 − ) (2.96)
𝑅 𝑅 𝑅

Buffer layer (region that assembles the viscous dominated layer near the
wall with the turbulent region), van Driest proposed (Schlichting, 1979):

1 𝑦𝑢∗
𝑙 = 𝑦 [1 − exp (− )] , 𝐴 = 26 (2.97)
𝐴 

:: 63 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

3.- VELOCITY DISTRIBUTIONS AND FRICTION FACTORS IN


TURBULENT FLOWS

Smooth or rough wall?


Before presenting the velocity distribution, it is necessary to define the
type of wall is the solid boundary of the conduit. Any surface observed with the
proper magnifier will look rough, as sketched in Fig. 3.1. In fluid mechanics, the
definition of smooth or rough wall does not depend only on the size of the

Fig. 3.1.- Roughness in a surface


roughness, but also in the thickness of the viscous sublayer. The roughness size
is commonly denoted by 𝜀 and the viscous sublayer thickness by 𝛿𝑉 . Thus, if the
roughness is large enough to avoid the presence of a viscous sublayer (𝜀 > 𝛿𝑉 )
the wall is hydraulically rough. On the contrary, if the viscous sublayer covers
all the roughness (𝜀 < 𝛿𝑉 ) the wall is hydraulically smooth. When the average
size of the roughenss is around the viscous sublayer thickness (𝜀 ~ 𝛿𝑉 ), the wall
is in transition smooth-rough. To determine the limits that define the type of
wall, we will follow Prandtl analysis to define the constant 𝐶 of the logarithmic
velocity distribution (Schlichting, 1979, p.589). As it is observed in Fig. 2.5, there
is a distance 𝑦0 at which the logarithmic profile is equal to cero (𝑢̅ = 0). Taken
it as a boundary condition, we obtain from Eq. 2.95:

𝑢∗
𝑢̅ = (ln(𝑦) − ln(𝑦0 )) (3.1)

The distance 𝑦0 is of the order of magnitude of 𝛿𝑉 . Using dimensional
analysis, we can get a dimensionless parameter involving 𝑦0 . The relevant
variables near the bottom is the fluid properties (viscosity 𝜇 and density 𝜌), the
shear stress acting on the wall (𝜏0 ). Thus, we can expect a functional relationship
of the form 𝑦0 = 𝑓(𝜇, 𝜌, 𝜏0 ). The number of variables is 𝑛 = 4 and the number of
fundamental dimensions is 𝑟 = 3 (F,L,T). Applying the Buckingham  theorem,
we have only one dimensionless parameter (𝑛 − 𝑟 = 1), given by:

𝜏01 2 𝜌1⁄2 𝑦0
= (3.2)
𝜇

:: 64 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Multiplying and dividing by 𝜌1⁄2 and recalling 𝑢∗ = √𝜏0 ⁄𝜌 and  = 𝜇 ⁄𝜌 we


get:

𝑢∗ 𝑦0
= (3.3)

Thus we have a functional relationship () = 0, that means that  must
be a constant. Calling 𝛽 such a constant:

𝑢∗ 𝑦0
=𝛽 (3.4)

Replacing 𝑦0 from Eq. 3.4 into Eq. 3.1 and dividing by the shear velocity:

𝑢̅ 1 𝑦𝑢∗
= (ln ( ) − ln(𝛽)) (3.5)
𝑢∗  
It is customary to define the inner (or wall) variables:

𝑢̅ 𝑦𝑢∗
𝑢̅+ ≡ , 𝑦+ ≡ (3.6)
𝑢∗ 
Eq. 3.5 is re-written as:

1
𝑢̅+ = (ln(𝑦 + ) − ln(𝛽)) (3.7)

The logarithmic profile written in dimensionless form as Eq. 3.7 will help
us to define the viscous sublayer and buffer region thicknesses. However, before
doing that, it will be useful to know the velocity distribution in the viscous
sublayer.
Velocity distribution in the viscous sublayer
The characteristic of the viscous sublayer is that the viscous stresses are
much larger than the turbulent ones, thus we can write 𝑇𝑦𝑥 = 𝜏𝑉𝑦𝑥 . In addition,
as this is a thin layer, we can assume that the shear stress is not so different
than its value on the bottom, 𝜏0 . Thus, we can write:

𝑑𝑢̅
𝜇 = 𝜏0 (3.8)
𝑑𝑦

Integrating and imposing the boundary condition 𝑦 = 0, 𝑢̅ = 0:

𝜇 𝑢̅ = 𝜏0 𝑦 (3.9)

:: 65 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Dividing Eq. 3.9 by 𝜌:

𝜇 𝜏0
𝑢̅ = 𝑦 (3.10)
𝜌 𝜌

In terms of the inner variables:

𝑢̅+ = 𝑦 + (3.11)

Thus, in the viscous sublayer the velocity varies linearly with the distance
to the wall.
The viscous sublayer thickness
At this point, we have done a lot to determine the viscous sublayer
thickness, but we need to rely on experimental data. The utility of having the
velocity distribution in terms of dimensionless variables is that it is possible to
compare data obtained under different flow conditions. Some experimental data
is presented in Fig. 3.2, which was taken from White (1991, Fig. 6.11). In
addition to the experimental data, the linear velocity profile form the viscous
sublayer (Eq. 3.11), and the logarithmic profile form the turbulent region
(Eq.3.5) are also plotted.
Obviously, the thickness of the viscous sublayer is defined by the distance
from the wall where the experimental data departs from the linear profile (that
appears as a curve in a semilogarithmic plot). This happens at 𝑦 + ≈ 5. Thus, the
dimensionless viscous sublayer thickness, 𝛿𝑉+ = 𝛿𝑉 𝑢∗ ⁄ is taken as:

𝛿𝑉+ = 5 (3.12)

The experimental data join the line defined by the logarithmic distribution
around 𝑦 + ≈ 70. That means that for 𝑦 + > 70 turbulence dominates over viscous
effects. Finally the buffer layer, or region where the viscous effects are as
important as the turbulent ones is comprised in 5 < 𝑦 + < 70. Often, the buffer
region is neglected and 𝛿𝑉+ ~ 11 is taken. Below that height the motion is
considered laminar and above, fully turbulent.
For the flow in a smooth conduit, the three regions mentioned before are
present. We are still talking about smooth walls but we have not given a precise
definition regarding when a wall can be considered smooth. Following Prandtl
(Schlichting, 1979) and von Karman (1930), a wall is hydraulically smooth when
the roughness is completely contained within the viscous sublayer. A wall is
hydraulically rough when the roughness size is large enough to preclude the

:: 66 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

̅+
𝒖

LINEAR DISTRIBUTION

LOGARITHMIC PROFILE

𝒚+
Fig. 3.2.- Experimental data and velocity distributions for the viscous
sublayer (Eq. 3.11) and turbulent region (Eq. 3.7).

existence of a region. As the roughness is much greater than the viscous


sublayer, the resistance comes mainly from the drag of the protruding elements
on the surface. The conceptual idea of both kind of walls is sketched in Fig. 3.3.
There is a transition from the hydraulically smooth to the rough surface, in
which both viscous stress and drag are important in the generation of flow
resistance. This originates a wall in the transition smooth-rough regime.
A question regarding how the roughness size is determined. It is not
measured directly. Actually it is an equivalent roughness. The size 𝜀 is the result
of a distribution of sizes that gives the same resistance to the flow than that
obtained by Nikuradse in experiments with uniform sand grains of size 𝑘𝑆 glued
to the inner wall of pipes, as shown in Fig. 3.4. Although both kind of roughness
have the same effects for the hydraulically smooth and rough walls, they present
a difference in the transition. For Nikuradse’s uniform grains, there is a sudden
passage from smooth to rough wall. In the case of 𝜀, there is a gentle transition
from smooth to rough due to the several sizes involved in the distribution
roughness heights. The value 𝜀 corresponds to that obtained for commercial
pipes.

:: 67 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝛿𝑉

a) SMOOTH WALL

𝛿𝑉 𝛿𝑉
𝛿𝑉
𝜀

b) ROUGH WALL
Fig. 3.3.- Definition of hydraulically smooth and rough surfaces.
a) On a smooth wall, the viscous sublayer cover completely the
roughness of the wall. The flow in the turbulent region is not
influenced by the size of the roughness and the flow resistance is
due to the viscous shear acting in the viscous sublayer.
b) The viscous sublayer is completely destroyed by the roughness
elements and the flow resistance is originated by the drag exerted
by them.

𝑘𝑆

Fig. 3.4.- Nikuradse’s roughness. Sand grains were glued to the wall of the
pipes.
Based on the limits of the different regions indicated before, the
hydrodynamic type of wall is defined depending on which of them the roughness
is contained. It is expressed in terms of the roughness size made dimensionless
with the inner variables:

𝑘𝑆 𝑢∗
𝑘𝑆+ = (3.13)

Due to the arbitrariness in the definition of the region boundaries, the
values assigned to 𝑘𝑆+ to define the type of wall vary slightly among different
authors, as shown in Table 3.1,although those assigned by Prandtl are the most
used in practice.

:: 68 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Table 3.1.- DIMENSIONLESS ROUGHNESS SIZE AND HYDRAULIC TYPE


OF WALL

PRANDTL WHITE
HYDRAULIC TYPE OF WALL
(Schlichting, 1979) (White, 1991)

SMOOTH 𝑘𝑆+ < 5 𝑘𝑆+ < 4

TRANSITION SMOOTH -ROUGH 5 ≤ 𝑘𝑆+ ≤ 70 4 ≤ 𝑘𝑆+ ≤ 60

ROUGH 𝑘𝑆+ > 70 𝑘𝑆+ > 60

VELOCITY DISTRIBUTION IN THE FULLY TURBULENT REGION


(𝒚+ > 𝟕𝟎)
As it was indicated, the velocity distribution in the fully turbulent region
is directly related to the flow resistance mechanism, which depends on the type
of wall. Thus, the distance 𝑦0 at which the boundary condition is evaluated
should consider it. Eq. 3.1 can be written as:

𝑢̅ 1 𝑦
= ln ( ) (3.14)
𝑢∗  𝑦0

The type of wall should be taken into account in the value of 𝑦0 .

Hydraulically smooth wall (𝒌+


𝑺 < 𝟓)

In this case, we already found that 𝑦0 ~ ⁄𝑢∗ (Eq. 3.4) and Eq. 3.7 is valid,
which usually is written as:

1
𝑢̅+ = ln(𝑦 + ) + 𝐵 (3.15)

 and 𝐵 are determined from experimental data, such as those presented
in Fig. 3.2, resulting  = 0.4 and 𝐵 = 5.5 (Nikuradse, 1933).

Hydraulically rough wall (𝒌+


𝑺 > 𝟕𝟎)

The flow resistance mechanism is due to the drag generated by the grains
large enough to preclude the existence of a viscous sublayer. In this case, it is
relevant the size of the roughness, so 𝑦0 ~ 𝑘𝑆 and Eq. 3.14 becomes:

𝑢̅ 1 𝑦
= ln ( ) + 𝐶 (3.16)
𝑢∗  𝑘𝑆

:: 69 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝐶 is determined from experimental measurements like that presented in


Fig. 3.5, which correspond to the data taken by Nikuradse (1933), and takes the
value 8.5.

Fig. 3.5.- Velocity distribution in a hydraulically rough pipe. Note that in


𝒚
the horizontal axis is 𝐥𝐨𝐠 𝟏𝟎 ( ) (Nikuradse, 1933).
𝒌𝑺

Wall in smooth-rough transition (𝟓 ≤ 𝒌+


𝑺 ≤ 𝟕𝟎)

I this case, the flow resistance is the result of both viscous and drag effects.
Thus, a functional relationship 𝑦0 = 𝑓(𝜇, 𝜌, 𝜏0 , 𝑘𝑆 ) should exist. From the
Buckingam  theorem we obtain:

𝑦0 𝑢∗ 𝑘𝑆
1 = , 2 = (3.17)
𝑘𝑆 
The dimensionless relation between both dimensionless parameters is
written as:

𝑦0 +
( ,𝑘 ) = 0 (3.18)
𝑘𝑆 𝑆

:: 70 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Which is the same that:

𝑦0
= 1 (𝑘𝑆+ ) (3.19)
𝑘𝑆

Thus, 𝑦0 is given by:

𝑦0 = 𝑘𝑆 1 (𝑘𝑆+ ) (3.20)

Replacing 𝑦0 from Eq. 3.20 into Eq. 3.14:

𝑢̅ 1 𝑦
= ln ( ) + 𝐴(𝑘𝑆+ ) (3.21)
𝑢∗  𝑘𝑆

The function 𝐴(𝑘𝑆+ ) is determined computing the deviation shown by the


experimental data corresponding to the smooth-rough transition with respect to
the velocity distribution for the rough wall (Eq. 3.16), resulting the graphic
relationship shown in Fig. 3.6 (Nikuradse, 1933).

𝐴(𝑘𝑆+ )

Fig. 3.6.- Function 𝑨(𝒌+


𝑺 ) of the velocity distribution in fully turbulent
region of a pipe with wall in the transition smooth-rough regime. Note
that in the horizontal axis is 𝐥𝐨𝐠 𝟏𝟎 (𝒌+
𝑺 ) (Nikuradse, 1933).

HYDRAULICALLY SMOOTH WALL (𝒌+


𝑺 < 𝟓). Velocity distribution in the
buffer layer (𝟓 ≤ 𝒚 ≤ 𝟕𝟎)
+

In hydraulically smooth walls, a buffer layer matching the viscous


sublayer with the fully turbulent region exists. Spalding (1961) gives two

:: 71 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

expressions for the velocity profile in the buffer layer. The difference between
them is the presence of a term raised to the fourth power. However, the equation
without this term adjust very well to the experimental data, as shown by White
(1991, Fig. 6.11). The equation is:

−𝐵
(𝑢+ )2 (𝑢+ )3
[𝑒 𝑢
+
+
𝑦 =𝑢 +𝑒 +
− 1 − 𝑢 −
+
− ] (3.22)
2 6

Note that Eq. 3.22 is not explicit in 𝑢+ , with  = 0.4 and 𝐵 = 5.5. The range
of validity of this equation extends further than the buffer layer, spanning 0 ≤
𝑦 + < 300.
The different velocity profiles arising in the turbulent flow are
summarized in Table 3.2.

Table 3.2.- SUMMARY OF VELOCITY DISTRIBUTIONS IN TURBULENT


FLOWS
TYPE OF WALL VELOCITY DISTRIBUTION

Viscous sublayer, 𝑦 + < 5


𝑢+ = 𝑦 +

HYDRAULICALLY SMOOTH Buffer layer, 5 ≤ 𝑦 + ≤ 70


𝑘𝑆+ < 5 Eq. 3.22

Fully turbulent region, 𝑦 + > 70


1
𝑢+ = ln(𝑦 + ) + 5.5

1 𝑦
TRANSITION SMOOTH-ROUGH 𝑢+ = ln ( ) + 𝐴(𝑘𝑆+ )
 𝑘𝑆
5 ≤ 𝑘𝑆+ ≤ 70
𝐴(𝑘𝑆+ ) from Fig. 3.6

HYDRAULICALLY ROUGH 1 𝑦
𝑢+ = ln ( ) + 8.5
𝑘𝑆+ > 70  𝑘𝑆

:: 72 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

FRICTION COEFFICIENT OF PIPES WITH TURBULENT FLOWS


The friction coefficient 𝑓, necessary to compute the frictional energy loss
(Eq. 47), can easily be computed from the velocity distributions presented in the
previous sections. Eq. 68 gives the relation between 𝑓 and the cross-sectional
mean velocity, 𝑈:

𝑢∗ 𝑓 (68)
=√
𝑈 8
The velocity distributions have the form 𝑢̅⁄𝑢∗ as a function of the distance
from the wall. It looks evident that, integrating the velocity distribution across
the flow section will allows to obtain 𝑈⁄𝑢∗ , hence the friction factor 𝑓. The
analysis that follows will consider cylindrical pipes because in this case exists an
extensive set of data that permits to confirm the validity of the expressions
obtained. The results can be generalized to other geometries.

Friction factor for turbulent flow with hydraulically smooth walls


The velocity distribution in the fully turbulent region is given by Eq. 3.15:

1
𝑢̅+ = ln(𝑦 + ) + 𝐵 (3.15)

The cross-sectional mean velocity is computed from

1
𝑈= ∫ 𝑢̅ 𝑑𝐴 (3.23)
𝐴 𝐴
For a pipe with inner diameter 𝐷, 𝐴 = 𝜋𝐷2 ⁄4. The distance from the wall
(𝑦) is related to the radial distance from the axis of the pipe (𝑟) trough 𝑦 + 𝑟 =
𝐷⁄2. The element of surface is 𝑑𝐴 = 2𝜋𝑟𝑑𝑟. In order to simplify the computations,
we will neglect the viscous sublayer and buffer regions, i.e., we will assume that
Eq. 3.15 is valid in all the flow domain. Thus, from Eq.s 3.15 and 3.23:
𝐷
4 2 𝑢∗ 𝐷 𝑢∗ (3.24)
𝑈= ∫ ( ln (( − 𝑟) ) + 𝑢∗ 𝐵) 2𝜋𝑟𝑑𝑟
𝜋𝐷 𝑜 
2 2 

𝐷
8 𝑢∗ 2 𝐷 𝑢∗ (3.25)
𝑈 = 2 ∫ (ln ( − 𝑟) + ln ( ) + 𝐵) 𝑟𝑑𝑟
𝐷  𝑜 2 

:: 73 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝐷 𝐷
8 𝑢∗ 2 𝐷 𝑢∗ 2
𝑈 = 2 [∫ ln ( − 𝑟) 𝑟𝑑𝑟 + (ln ( ) + 𝐵) ∫ 𝑟𝑑𝑟] (3.26)
𝐷  𝑜 2  𝑜

𝐷
8 𝑢∗ 2 𝐷 𝑢∗ 1 𝐷 2
𝑈 = 2 [∫ ln ( − 𝑟) 𝑟𝑑𝑟 + (ln ( ) + 𝐵) ( ) ] (3.27)
𝐷  𝑜 2  2 2

A comment regarding the integral in Eq. 3.27. It will give a term


containing (𝐷 − 2𝑟)ln(𝐷 − 2𝑟)that when evaluated at 𝑟 = 𝐷⁄2 gives 0 × (−∞). To
overcome this problem, l’Hôpital rule must be applied. Grouping the sum of
logarithms in the logarithm of the product, the following result is obtained:

𝑢∗ 1 𝐷𝑢∗ 3
𝑈= [(ln ( ) − ) + 𝐵] (3.28)
 2  2

We can rearrange Eq. 3.28 to give the structure of Eq. 68:

𝑢∗ 𝑓 1
=√ = (3.29)
𝑈 8 1 (ln (1 𝐷𝑢∗ ) − 3) + 𝐵
 2  2
Generally, the flow mean velocity 𝑈 = 𝑄/𝐴 is known, and not the shear
velocity, and it is convenient to change the argument of the logarithm in Eq. 3.29
multiplying and diving it by 𝑈:

𝐷𝑢∗ 𝐷𝑈 𝑢∗ 𝑓 (3.30)
= = 𝑅𝑒√
  𝑈 8

Replacing Eq. 3.30 in Eq. 3.29 and working with its reciprocal:
1 1 3
= [ln(𝑅𝑒√𝑓) − ln(√32) − + 𝐵] (3.31)
√𝑓 √8 2

Using the values  = 0.4 and 𝐵 = 5.5:


1
= 0.884ln(𝑅𝑒√𝑓) − 0.913 (3.32)
√𝑓
Eq. 3.32 is usually expressed in terms of the decimal logarithm:

:: 74 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

1
= 2.035log(𝑅𝑒√𝑓) − 0.913 (3.33)
√𝑓
Eq. 3.33 was obtained by Prandtl in 1935 (White, 1991). As the buffer
region and viscous sublayer were neglected in its deduction, the numerical
constants that appear in it were re-computed in order to fit the experimental
data taken by Nikuradse (1933), resulting:
1
= 2log(𝑅𝑒√𝑓) − 0.8 (3.34)
√𝑓
This is Prandtl’s universal law of friction for smooth pipes (Schlichting,
1979) and it is not have limitations regarding the Reynolds number (within the
turbulent regime). Computation of 𝑓 from Eq. 3.34 must be done graphically or
numerically, usually by iteration. A simpler equation is the empirical
relationship presented by Blasius in 1913, but limited to 3 × 103 ≤ 𝑅𝑒 ≤ 105
(Blasius, 1913, p. 12)

0.3164
𝑓= (3.35)
𝑅𝑒 1⁄4
Friction factor for turbulent flow with hydraulically rough walls
The friction factor for turbulent flows with rough walls is obtained from
the velocity distribution given by Eq. 3.16.

𝑢̅ 1 𝑦
= ln ( ) + 𝐶 (3.16)
𝑢∗  𝑘𝑆

Assuming eq. 3.16 valid in all the cross-section, computing the average
velocity 𝑈 , rearranging the final result to form √𝑓 = √8 𝑢∗ ⁄𝑈, and using  = 0.4
and 𝐶 = 8.5, the following expression is obtained:
1 𝐷
= 0.942 ln (
2𝑘𝑆
) + 1.68 (3.36)
√𝑓
In terms of the logarithm in base 10:
1 𝐷
= 2.2 log ( ) + 1.68 (3.36)
√𝑓 2𝑘𝑆

The numerical constants are adjusted to the experimental data, resulting:


1 𝐷
= 2 log (
2𝑘𝑆
) + 1.74 (3.37)
√𝑓

:: 75 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Frequently, the additive constant 1.74 is introduced into the logarithm


argument, giving:
1 𝐷
= 2 log (3.7
𝑘𝑆
) (3.38)
√𝑓

Friction factor for turbulent flow with walls in transition smoot-rough


An analytical derivation of a relationship for the friction factor for the
turbulent flow with walls in the transition smooth-rough regime is not possible
because we do not have a simple analytical equation for 𝑢̅⁄𝑢∗ . In this case, the
velocity profile is given by Eq. 3.21, with the additive constant depending of 𝑘𝑆+ ,
given in Fig. 3.6. In order to determine a relation for 𝑓, we will analyse the
deviation of the experimental data from the friction factor associated to a
hydraulically rough wall. In order to avoid confusion with the friction factors, we
will use the sub-indexes 𝑆 to denote the smooth wall and 𝑅 for the rough. We can
write for smooth wall friction coefficient (Eq. 3.34):
1
2log(𝑅𝑒√𝑓𝑆 ) − = 0.8 (3.39)
√𝑓𝑆

𝑈𝐷 𝑢∗ 1
2log (

√8
𝑈
)− = 0.8 (3.40)
√𝑓𝑆

𝑢∗ 𝐷 1
2log ( √8) − = 0.8 (3.41)
 √𝑓𝑆

𝑢∗ 𝑘𝑆 𝐷 1
2log ( √8) − = 0.8 (3.42)
 𝑘𝑆 √𝑓𝑆

𝐷 1
2log ( ) + 2log(𝑘𝑆+ ) + 2log(√8) − = 0.8 (3.43)
𝑘𝑆 √𝑓𝑆

:: 76 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝐷 1
2log (3.7 ) − 2log(3.7) + 2log(𝑘𝑆+ ) + 2log(√8) − = 0.8 (3.44)
𝑘𝑆 √𝑓𝑆

We identify the first term of Eq. 3.44 as 1⁄√𝑓𝑅 (Eq. 3.38):

1 1
− 2log(3.7) + 2log(𝑘𝑆+ ) + 2log(√8) − = 0.8 (3.45)
√𝑓𝑅 √𝑓𝑆

1 1
− = 0.8 + 2log(3.7) − 2log(√8) − 2log(𝑘𝑆+ ) (3.46)
√𝑓𝑅 √𝑓𝑆

1 1
− = 1.033 − 2log(𝑘𝑆+ ) (3.47)
√𝑓𝑅 √𝑓𝑆

Eq. 3.47 is the deviation (in terms of 1⁄√𝑓 ) of the friction factor for smooth
walls from the friction factor for rough walls. Now, we will defining the deviation
of the friction factor for any type of wall from that associated to the rough surface
as:
𝐷 1
𝑀 = 2log (3.7
𝑘𝑆
)− (3.48)
√𝑓
As Eq. 3.48 is valid for any type of wall, 𝑀 should be function of 𝑘𝑆+ . Using
the experimental data presented by Nikuradse (1933) is possible to know 𝑀(𝑘𝑆+ ).
In effect, with the data is possible to generate the Table. 3.3. The first four
columns contains Nikuradse’s data and the fifth is computed using Eq. 3.48.

Table 3.3.- ORGANIZATION OF NIKURADSE’S DATA TO GENERATE THE


FUNCTION 𝑀(𝑘𝑆+ )
𝐷 𝑘𝑆 𝑢∗
𝑅𝑒 𝑓 𝑘𝑆+ = 𝑀
𝑘𝑆 
: : : : :
: : : : :
: : : : :

:: 77 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

The relation between 𝑀 and 𝑘𝑆+ is presented in Fig. 3.7, taken from the
article by Colebrook and White (1937) and corresponds to line labelled
“Nikuradse sanded”. The straight line labelled “smooth law” is Eq. 3.47.
Obviously, for rough walls 𝑀 = 0 (“rough law” in the figure).

𝑘𝑆+
Fig. 3.7.- Deviation of 𝟏⁄√𝒇 from the rough wall. The function 𝑴
covers all type of walls.

The objective of the article by Colebrook and White (1937) was “to
determine how the nature of the roughness influenced the transition”. To
accomplish that goal, sand of different sizes and following some particular
arrangements were glued to the pipe. But they also added information of
commercial pipes reported by other authors (Freeman, quoted by Mills, 1923;
and Heywood, 1924). The commercial pipes generate the curve labelled as
“Galvanised and new wrought iron” in Fig. 3.7. The different behaviour among
the measurements carried out in commercial pipes and those obtained by
Nikuradse is result of the non-uniformity of the roughness size in the commercial
pipes. In this case, protuberances larger than the average disturb and destroy
the viscous sublayer before that the rest of the protuberances. This is a gradual
process that makes an earlier and gradual separation of the 𝑀 curve from the
smooth case when compared with pipes of uniform roughness.
As it was indicated before, it is customary to denote by 𝜀 the equivalent
roughness that arise in pipes with non-uniform size of protuberances. A
dimensionless equivalent roughness is defines as:

𝜀𝑢∗
𝜀+ ≡ (3.49)

:: 78 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Thus, for commercial pipes, the experimental relationship between 𝑀 and


𝜀 is found to be:
+

3.29
𝑀 = 2 log (1 + ) (3.50)
𝜀+
It is easy to verify that Eq. 3.50 does not have limitations regarding the
type of wall, covering all of them. For smooth walls, 𝜀 + is very small, and
3.29⁄𝜀 + ≫ 1, and 𝑀 ≈ 2 log(3.29⁄𝜀 + ) = 2 log(3.29) − 2 log(𝜀 + ), that is to say:

𝑀 = 1.033 − 2log(𝜀 + ) (3.51)

Eq. 3.51 is in agreement with Eq. 3.47, the deviation for smooth surfaces.
In the same way, for hydraulically rough walls, 𝜀 + is large, resulting 3.29⁄𝜀 + ≪
1, and 𝑀 ≈ 2log(1) = 0, i.e., there is not deviation from the friction factor for
rough surface.
Thus, it is possible to determine an expression for the friction factor for
turbulent flows, valid for all type of walls. In effect, replacing Eq. 3.50 into eq.
3.48:
3.29 𝐷 1
2 log (1 +
𝜀+
) = 2 log (3.7
𝜀
) − (3.52)
√𝑓
But

𝜀𝑢∗ 𝜀 𝐷𝑢∗ 𝜀 𝐷𝑈 𝑢∗ 𝜀 𝑓 (3.53)


𝜀+ = = = = 𝑅𝑒√
 𝐷  𝐷  𝑈 𝐷 8
Replacing Eq. 3.53 in Eq. 3.52 results:

1 3.29√8 𝐷 𝐷
− = 2 log (1 + ) − 2 log (3.7 ) (3.54)
√𝑓 𝑅𝑒√𝑓 𝜀 𝜀

1 3.29√8 𝐷 𝜀
− = 2 log [(1 + ) ] (3.55)
√𝑓 𝑅𝑒√𝑓 𝜀 3.7𝐷

1 𝜀 3.29√8 1
− = 2 log ( + ) (3.56)
√𝑓 3.7𝐷 3.7 𝑅𝑒√𝑓
Finally, the friction factor for turbulent flows is pipes is given by:

:: 79 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

1 𝜀 2.51
= −2 log ( + ) (3.57)
√𝑓 3.7𝐷 𝑅𝑒√𝑓

Eq. 3.57 was presented in an article written by Colebrook in 1939.


Although he is the only author of the article, Colebrook recognizes White's
contribution in the development of the formula, which is known as the
Colebrook-White equation. Eq. 3.57 is not explicit in 𝑓, therefore computation
must be done iteratively, which is not a problem with the calculators existing
nowadays. However, in the 1940s, it resulted a cumbersome task when it had to
be computed manually o with slide rules. Thus, the graphic relationship shown
in Fig. 3.8 was presented by Rouse (1943). The graph contains a set of curves
with the relative roughness (𝜀 ⁄𝐷) as parameter. Depending on the known
variables, the diagram can be accessed from the lower horizontal axis, which
contains 𝑅𝑒√𝑓. Considering the Darcy-Weisbach equation (Eq.47), 𝑓 = 2𝐽𝑔𝐷 ⁄𝑈 2 ,
the term 𝑅𝑒√𝑓 can be written as: 𝑅𝑒√𝑓 = 𝑈𝐷⁄ √2𝐽𝑔𝐷 ⁄𝑈 2 = √2𝐽𝑔 𝐷 3⁄2 ⁄, which
is independent of the velocity. Thus, if the viscosity, gradient of the energy loss
( 𝐽), pipe diameter and roughness are known the velocity (and hence the
discharge) are determined directly. On the other hand, if the velocity and pipe
characteristics are known and the friction factor has to be computed, the
Reynolds number can be determined and the graph is accessed from the upper
horizontal axis until reach the line corresponding to the relative roughness and
𝑓 (or 1⁄√𝑓 ) is read in the vertical axis. Note that 𝑅𝑒 follow curved lines. (In
Rouse’s notation, 𝑆 = 𝐽).
However, the most known graphic relationship is that due to Moody
(1944), who only differs from Rouse’s in the principal axis. Moody used the
Reynolds number in the lower horizontal axis. Thus, we have to follow a vertical
straight line for a given 𝑅𝑒 (and not curved lines as in Rouse diagram). Note that
Moody indicates that the units of velocity, diameter and kinematic viscosity
must be ft/s, ft and ft2/s, respectively. Obviously, this is not necessary. As 𝑅𝑒 is
dimensionless, the only requirement is that the units must be in any coherent
system of units.
Note that we have presented relationships for the friction factor for two
flow regimes: laminar (𝑓 = 64⁄𝑅𝑒 ) and turbulent (Eq. 3.57). We have not
presented a relation for the transition laminar-turbulent regime because it is
very uncertain and small variations of 𝑅𝑒 mean large variations in 𝑓. In general,
flows in the transition regime are not usual in civil or environmental
engineering.

:: 80 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Fig. 3.8.- Rouse diagram (Rouse, 1943)

A large amount of equations for 𝑓 in the turbulent regime can be found in


the literature. Many of them are experimental relationships determined for
specific pipes and materials. Also, there a many relationships that are fittings to
the Moody diagram (Eq.3.57), with only goal to have explicit expressions for 𝑓.
Some of them are the following (Beluco and Camano, 2016):
Haaland:
1 𝜀 1.1 6.9
= −1.8 log [( ) + ] (3.58)
√𝑓 3.7𝐷 𝑅𝑒

Barr:
1 𝜀 5.15
= −2 log [ + 0.892 ] (3.59)
√𝑓 3.7𝐷 𝑅𝑒

Eck:
1 𝜀 15
= −2 log [ + ] (3.60)
√𝑓 3.715𝐷 𝑅𝑒

:: 81 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

As the above equations (an others similar) are approximations to Eq. 3.57,
the friction factor computed with them will have an error.
Finally in Table 3.4, typical values of the roughness size for some
materials are presented.
Table 3.4.- TYPICAL VALUES OF THE ROUGHNESS SIZE
MATERIAL 𝜀 (mm)

Riveted steel 0.92 - 9.2


Concrete 0.31 - 3.1
Ductile iron 2.6
Wood stave 0.09 - 0.18
Galvanized iron 0,15
Cast iron – asphalt dipped 0,12
Cast iron uncoated 0,25
Carbon steel or wrought iron 0,045
Stainless steel 0,045
Fiberglass 0,005
Drawn tubing – glass, brass, plastic 0,0015
Copper 0,0015
Aluminium 0,0015
PVC 0,0015
Red brass 0,0015

:: 82 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Fig. 3.9.- Moody’s digram (Moody, 1944)

:: 83 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

RESISTANCE LAWS IN OPEN CHANNEL FLOWS


The application of relationships like Eq. 3.57 has been extended to the
flow in open channels, adjusting the numerical coefficients with ad hoc
experimental data (Yen, 2002), and replacing the pipe diameter 𝐷 by the
hydraulics radius 𝑅𝐻 = 𝐴⁄ , where 𝐴 is the flow area and  is the wetted
perimeter. Eq. 3.57 is transformed into:

1 𝑘𝑆 𝐾3
= −𝐾1 log ( + ) (3.61)
√𝑓 𝐾2 𝑅 4𝑅𝑒𝑅𝐻 √𝑓

In Eq. 3.61, 𝑅𝑒𝑅𝐻 is the Reynolds number based on the hydraulic radius, i.e.,
𝑅𝑒𝑅𝐻 = 𝑈𝑅𝐻 ⁄ . 𝐾1 , 𝐾2 and 𝐾3 are coefficients fitted by different authors. Yen
(2002) presents the summary shown in Table 3.5.

Table 3.5.- COEFFICIENTS OF EQ. 3.57 ACCORDING TO DIFFERENT


AUTHORS (YEN, 2002)

CHANNEL
REFERENCE 𝐾1 𝐾2 𝐾3 REMARKS
GEOMETRY
Full circular pipe Colebrook (1939) 2.0 14.83 2.52

Wide channel Keulegan (1938) 2.03 11.09 3.41

Wide channel Rouse (1946, p. 214) 2.03 10.95 1.70

Wide channel Thijsse (1949) 2.03 12.2 3.033


Sayre and Albertson
Wide channel 2.14 8.888 7.17
(1961)
Wide channel Henderson (1966) 2.0 12.0 2.5

Wide channel Graf (1971, p. 305) 2.0 12.9 2.77

Wide channel Reinius (1961) 2.0 12.4 3.4


Width/depth
Rectangular Reinius (1961) 2.0 14.4 2.9
=4
Width/depth
Rectangular Reinius (1961) 2.0 14.8 2.8
=2
Rectangular Zegzhda (1938) 2.0 11.55 0 Dense sand

:: 84 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

Table 3.6.- EQUIVALENT SIZE ROUGHNESS ACCORDING TO


DIFFERENT AUTHORS (YEN, 2002)

RESEARCHER CHARACTERISTIC SIZE 𝛼𝑠


Ackers and White (1973) 𝑑35 1.23
Strickler (1923) 𝑑50 3.3
Keulegan (1938) 𝑑50 1
Meyer-Peter and Muller (1948) 𝑑50 1
Thompson and Campbell (1979) 𝑑50 2.0
Hammond et al. (1984) 𝑑50 6.6
Einstein and Barbarossa (1952) 𝑑65 1
Irmay (1949) 𝑑65 1.5
Engelund and Hansen (1967) 𝑑65 2.0
Lane and Carlson (1953) 𝑑75 3.2
Gladki (1979) 𝑑80 2.5
Leopold et al. (1964) 𝑑84 3.9
Limerinos (1970) 𝑑84 2.8
Mahmood (1971) 𝑑84 5.1
Hey (1979), Bray (1979) 𝑑84 3.5
Ikeda (1983) 𝑑84 1.5
Colosimo et al. (1986) 𝑑84 3-6
Whiting and Dietrich (1990) 𝑑84 2.95
Simons and Richardson (1966) 𝑑85 1
Kamphuis (1974) 𝑑90 2.0
van Rijn (1982) 𝑑90 3.0

When Eq. 3.61 is applied to natural channels (with fixed boundaries), it is also
necessary to define the roughness, 𝑘𝑆 . The bed of natural channels, when it is
formed by granular non-cohesive sediments, usually contains a wide range of
particle sizes. The matter of what is the appropriate size to represent 𝑘𝑆 depends
on the researcher, and usually is expressed in terms of the diameter 𝑑𝑋 of the
size distribution (𝑑𝑋 is the size of the particle under which it is found the 𝑋% of
the sediment):

:: 85 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

𝑘𝑆 = 𝛼𝑆 𝑑𝑋 (3.62)

Yen (2002) gives a list of the diameter 𝑑𝑋 used by different authors and
the corresponding value of the coefficient 𝛼𝑆 , which is reproduced in Table 3.6.

:: 86 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

REFERENCES
BELUCO, A. and E.B. CAMANO SCHETTINI (2016) “An Improved Expression
for a Classical Type of Explicit Approximation of the Colebrook White Equation
with Only One Internal Iteration”, International Journal of Hydraulic
Engineering, Vol. 5, No. 1, pp. 19-23.
BLASIUS, H. (1913) “Das Aehnlichkeitsgesetz bei Reibungsvorgangen in
Flussigkeiten”, Mitteilungen Über Forschungsarbeiten auf dem Gebiete des
Ingenieurwesens, Heft 131, pp. 1-40.
BOUSSINESQ, J. (1877) Essai sur la théorie des eaux courantes, Mémoires
présentés par divers savants à l’Académie des Sciences. Imprimerie Nationale,
Paris.
BRADSHAW, P. (1972) “The understanding and prediction of turbulent flow”,
Aeronautical Journal, Vol. 76, No. 739, pp. 403-418.
COLEBROOK, C.F. (1939) “Turbulent Flow in Pipes, with particular reference
to the Transition Region between the Smooth and Rough Pipe Laws”, J. Inst.
Civil Eng. (London), Vol. 11, No. 4, pp. 133–156.
COLEBROOK, C.F. and C.M. WHITE (1937) “Experiments with Fluid Friction
in Roughened Pipes”, Proc. R. Soc. Lond. A, Vol. 161, pp. 367-381.
CHOW, V.T. (1988) Open-channel Hydraulics, McGraw-Hill.
HINZE, J.O. (1975) Turbulence, 2nd edition, McGraw‐Hill, New York.
JACKSON, D., and B. LAUNDER (2007) “Osborne Reynolds and the Publication
of His Papers on Turbulent Flow” Annu. Rev. Fluid Mech. Vol. 39, pp.19–35.
MOODY, L.F. (1944). "Friction factors for pipe flow." Trans. ASME, Vol. 66, pp.
671-678
MOSKVITCH, K. (2014) “Fiendish million-dollar proof eludes mathematicians.
Proposed solution to Navier–Stokes equations, used to model fluids, is wrong”,
Nature News (05 August 2014) doi:10.1038/nature.2014.15659
NAVIER, C.L.M.H. (1823) "Memoire sur les lois du mouvement des fluides,"
Mem. Acad. R. Sci. Paris, Vol. 6, pp. 389-416.
NIKURADSE, J. (1933) "Stromungsgesetze in rauhen Rohren." VDI-
Forschungsheft 361. Beilage zu "Forschung auf dem Gebiete des
Ingenieurwesens" Ausgabe B Band 4, July/August 1933. English translation:
“Laws of flow in rough pipes”, NACA TM 1292, Washington, November 1950.

:: 87 ::
Università degli Studi della Campania Luigi Vanvitelli
Dipartimento di Ingegneria Civile Design Edilizia e Ambientale

PRANDTL, L. (1925) “Bericht ber Untersuchungen zur ausgebildeten”


Turbulenz, Z. Angew. Math. Mech., Vol. 5, pp. 136-139.
REYNOLDS, O. (1883) “An experimental investigation of the circumstances
which determine whether the motion of water shall be direct or sinuous, and of
the law of resistance in parallel channels”. The Philosophical Transactions of the
Royal Society, Vol. 174, pp. 935-982.
REYNOLDS, O. (1895) “On the Dynamical Theory of Incompressible Viscous
Fluids and the Determination of the Criterion”, Philosophical Transactions of
the Royal Society of London. A, Vol. 186. (1895), pp.123-164.
ROUSE, H. (1943). "Evaluation of boundary roughness." Proceedings Second
Hydraulics Conference, Univ. of Iowa Studies in Engrg., Bulletin No. 27, pp. 105-
116.
SCHLICHTING, H. (1979) Boundary Layer Theory. 7 ed., McGraw-Hill, New
York.
SPALDING, D.B. (1961) “A Single Formula for the N"Law of the Wall””, J. appl.
Mech., Vol. 28, No. 3, pp. 455-457.
STOKES G.G. (1845) "On the Theories of Internal Friction of Fluids in Motion,"
Trans. Cambridge Phil. Soc., vol. 8, pp. 287-305. Also reprinted in Mathematical
and Physical Papers by G.G. Stokes, Vol. I., (1880), Cambridge at The University
Press. pp. 75-105.
TENNEKES, H. and J.J. LUMLEY (1972) A first course in turbulence. The MIT
Press, Massachusetts.
von KARMAN; T. (1930) “Mechanisclhe Ähnlchkeit und Turbulenz”, Nach-
richten von der Gesellschaft der Wisserschaften zu Göttingen, Fchgruppe I
(Mathematik) No. 5, pp. 58-76. English translation “Mechanical Similitude and
Turbulence”, NACA TM 611, Washington, March, 1931.
von KARMAN, T. (1937). “Turbulence”. The Journal of the Royal Aeronautical
Society, 41(324), 1109-1143. doi:10.1017/S0368393100103785
WHITE, F. (1991) Viscous Fluid Flow, 2nd Edition, McGraw-Hill, Inc., Singapore.
YEN, B.C. (2002) “Open Channel Flow Resistance”, J. Hyd. Engng., ASCE, Vol.
128, No. 1, pp. 29-39

:: 88 ::

You might also like