**Teaching and Learning of Fluid Mechanics**

Special Issue Editor **Ashwin Vaidya**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Special Issue Editor* Ashwin Vaidya Montclair State University USA

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Fluids* (ISSN 2311-5521) (available at: https://www.mdpi.com/journal/fluids/special issues/teaching learning fluids).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03936-443-5 (Pbk) ISBN 978-3-03936-444-2 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**



## **About the Special Issue Editor**

**Ashwin Vaidya** is a faculty member in the Department of Mathematics and Department of Physics & Astronomy at Montclair State University. He also directs the Complex fluids laboratory where students and faculty members work on experimental and theoretical problems pertaining to fluid mechanics. His general scholarly interests, besides fluid mechanics include complex systems, thermodynamics, math and science education and philosophy of science.

## *Editorial* **Teaching and Learning of Fluid Mechanics**

#### **Ashwin Vaidya**

Department of Mathematics, Montclair State University, Montclair, NJ 07043, USA; vaidyaa@mail.montclair.edu Received: 8 April 2020; Accepted: 10 April 2020; Published: 13 April 2020

Fluid mechanics is arguably one of the oldest branches of physics, and the literature on this subject is vast and complex. However, this subject has not sufficiently captured the interest of STEM educators like in other subjects such as quantum mechanics [1,2]. The objective of this collection, while not necessarily intended to generate education research, aims at bringing together various ways of teaching and learning about different topics in fluid mechanics.

Fluid mechanics occupies a privileged position in the sciences; it is taught in various science departments including physics, mathematics, environmental sciences and mechanical, chemical and civil engineering, with each highlighting a different aspect or interpretation of the foundation and applications of fluids. While scholarship in fluid mechanics is vast, expanding into the areas of experimental, theoretical and computational, there is little discussion among scientists about the different possible ways of teaching this subject or wide awareness of the how fluid mechanics plays a role in different disciplines. We believe there is much to be learned from an interdisciplinary dialogue about fluids for teachers and students alike.

The terms 'interdisciplinary', 'multidisciplinary' and 'transdisciplinary' have become common parlance in academia, but have been misunderstood and used without distinction [3]. Multidisciplinary teaching refers to diverse parallel viewpoints, with different goals and objectives being presented in the same setting while interdisciplinary or transdiscplinary refer to instances where goals overlap or unify completely [4]. In the context of education, inter- or transdisciplinary instruction allows students get to see the commonalities between different disciplines, thereby allowing students to make new meanings out of old ideas [5]. The education theorist, William Doll [5], articulates this idea very well: "Order emerges from interactions having just the 'right amount' of tension or difference or imbalance among the elements interacting." In a recent paper on education, my co-authors and I have argued that the synergy between different disciplines can result in the emergence of order, which we argue is nothing but creativity [6]. Doll's fluid analogy [5] for this idea is especially relevant to this issue:

"Emergence of creativity from complex flow of knowledge—example of Benard convection pattern as an analogy—dissipation or dispersal of knowledge (complex knowledge) results in emergent structures, i.e., creativity which in the context of education should be thought of as a unique way to arrange information so as to make new meaning of old ideas."

With this philosophy in mind, we have included all kinds of articles in this issue, including research on the pedagogical aspects of fluid mechanics, case studies or lesson plans at the undergraduate or graduate levels, articles on historical aspects of fluids, and novel and interesting experiments or theoretical calculations that can convey complex ideas in creative ways. The current volume includes 14 papers and showcases the work of scientists from different disciplines ranging from mathematics and physics to mechanical, environmental and chemical engineering. It truly is a wonderful collection and provides ideas on theoretical computational and experimental aspects of fluid mechanics that be implemented in a course in fluids in any department. The suitability of these papers ranges from early undergraduate to graduate level.

Overall, this issue contains papers in various, somewhat distinct categories. The articles [7–11] add and *reconsider fundamental ideas in fluid mechanics*. The paper by Brki´c and Praks [7] is devoted to the solution of Colebrook's friction equation. This basic equation is used to introduce interesting and sophisticated mathematical tools, such as the fixed point method and Padé approximation, among others. The article by Kariotoglou and Psillos [8] is a synthesis of previous research upon high school and undergraduate students on the teaching of concepts such as pressure in fluids. The two papers by Pal [9,10] are focused on fundamental ideas in thermomechanics; the papers discuss the ways in which important thermodynamic ideas can be introduced and elucidated in a fluid dynamics course through examples such as flow in pipes and flow through packed beds. Vianna et al. [11] take up the concept of *permeability* in porous media flow and discusses new computational concepts that can help convey such complex concepts.

Several collections in this issue deal with the *development of computational tools* to resolve important problems in fluid mechanics [12–14]. Addair and Jaeger [12] consider efficient and effective strategies to convey fundamental concepts in Computational Fluids Dynamics (CFD) to undergraduate students, such as the implementation of the finite volume method to solve Navier–Stokes equations. The paper by Battista [13] provides an open repository of several useful two-dimensional solvers written in MATLAB and Python 3. The contribution by Pawar and San [14] is about CFD Julia, a programming module that teaches the foundations of computational fluid dynamics (CFD). This piece is written for an upper-level undergraduate and early undergraduate course in fluid mechanics and uses the inviscid Burger's equation and the two-dimensional Poisson equation as examples.

Articles [10,15] treat *classical topics in fluid mechanics* in a very thorough and innovative manner and serve to guide the development of lectures on these topics in undergraduate and graduate courses. Medved, Davis and Vasquez [15] deal with the classical problem of a particle motion in a fluid. Pal's contribution [10] cited earlier regarding melding fluid mechanics and thermodynamics in the context of pipeline flow is yet another example of such a contribution.

Articles [16,17] provide *novel and creative new ways of introducing fluid mechanics* to students. They make rich connections between several disciplines and are guaranteed to capture the imagination of students. The paper by Mayer [16] discusses the seemingly simple example of an experiment on the flow of a fluid through a bottle. This seemingly mundane topic is a classic example of the richness and complexity of fluid mechanics, but also its ubiquity. Nita and Ramanathan [17] make a creative connection between fluids and music in their discussion of the physics of the Pan's flute. These are both truly exciting examples that can be introduced in undergraduate or graduate courses on fluids.

Articles [18–20] are more *pedagogically focused*. Huilier's contribution [18] is unique and provides a rich personal history of teaching fluid mechanics for over forty years at the University of Strasbourg. While the subject of fluid mechanics has evolved much in the last several decades, this reflective piece is a guide to all young instructors about how to adapt and evolve one's teaching and pedagogy. The paper by Potter and Wolff [19] discusses the experience and observations of interventions in a second-year fluid mechanics module, and provides insights into the teaching of fluid mechanics at their institution. A second goal of this paper is to conceptualize the use of Legitimation Code Theory (LCT) dimensions towards teaching strategies intended to facilitate improved learning outcomes. The LCT provides educators with a very useful tool to guide their curriculum planning and pedagogy. The article [20] by Valyrakis et al. focuses on the importance of lab-based fluid mechanics instruction. The authors demonstrate the value of using physical models ("floodopoly") to demonstrate complex geophysical processes such as sediment transport and flooding. A survey instrument is used to assess and demonstrate student understanding and perception of the subject.

#### **References**


© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Teaching Fluid Mechanics and Thermodynamics Simultaneously through Pipeline Flow Experiments**

#### **Rajinder Pal**

Department of Chemical Engineering, University of Waterloo, Waterloo, ON N2L 3G1, Canada; rpal@uwaterloo.ca; Tel.: +1-519-888-4567 (ext. 32985)

Received: 12 May 2019; Accepted: 28 May 2019; Published: 1 June 2019

**Abstract:** Entropy and entropy generation are abstract and illusive concepts for undergraduate students. In general, students find it difficult to visualize entropy generation in real (irreversible) processes, especially at a mechanistic level. Fluid mechanics laboratory can assist students in making the concepts of entropy and entropy generation more tangible. In flow of real fluids, dissipation of mechanical energy takes place due to friction in fluids. The dissipation of mechanical energy in pipeline flow is reflected in loss of pressure of fluid. The degradation of high quality mechanical energy into low quality frictional heat (internal energy) is simultaneously reflected in the generation of entropy. Thus, experiments involving measurements of pressure gradient as a function of flow rate in pipes offer an opportunity for students to visualize and quantify entropy generation in real processes. In this article, the background in fluid mechanics and thermodynamics relevant to the concepts of mechanical energy dissipation, entropy and entropy generation are reviewed briefly. The link between entropy generation and mechanical energy dissipation in pipe flow experiments is demonstrated both theoretically and experimentally. The rate of entropy generation in pipeline flow of Newtonian fluids is quantified through measurements of pressure gradient as a function of flow rate for a number of test fluids. The factors affecting the rate of entropy generation in pipeline flows are discussed.

**Keywords:** undergraduate education; fluid mechanics; pipeline flow; non-equilibrium thermodynamics; entropy generation; pressure loss; experimental studies

#### **1. Introduction**

Fluid mechanics, that is, the study of motion of fluids and forces in fluids, is relatively a less abstract subject as compared with thermodynamics. Students find it relatively easy to visualize and understand the motion of fluids and forces in fluids. For example, it is not difficult to convey to the students fluid mechanics concepts such as pressure, pressure distribution, viscosity, velocity distribution, velocity gradient, shear and normal stresses, mechanical energy dissipation and pressure loss in flow of fluids due to friction, etc. Most fluid mechanics quantities are directly measurable. For example, instruments are available to directly measure pressure or pressure drop, flow rate, local velocity, shear and normal stresses, etc. Thermodynamics, on the other hand, is a very abstract subject. Concepts such as entropy and entropy generation in real processes are illusive for students. There are no instruments available which can be used to directly measure entropy and entropy generation in real processes.

It is a well-known fact that students learn and understand concepts better and with relative ease through experiential learning. It is appropriate to quote here a famous Chinese saying credited to the Chinese philosopher Confucius, "I hear and I forget; I see and I remember; I do and I understand". To that end, the undergraduate laboratory experiments can play a very important role in providing experiential learning of concepts and theory to students.

In most undergraduate engineering programs (chemical, mechanical, civil, etc.) students are taught fluid mechanics through in-class instruction and laboratory experiments. In a typical undergraduate fluid mechanics laboratory, the students are required to do pipeline flow experiments involving measurement of pressure loss as a function of flow rate for different diameter pipes. The flow rate is measured with the help of a flowmeter and the pressure drop is measured with the help of pressure transducers. The test fluid is usually water and the experiments are carried out at room temperature. From pressure drop vs. flow rate experimental data, friction factor vs. Reynolds number data are calculated and compared with the available theoretical and empirical relations. What is completely missing in the undergraduate fluid mechanics experiments is the link between the pipeline flow experiments and the second law of thermodynamics, that is, entropy and entropy generation in pipeline flows. In order to appreciate and understand the second law of thermodynamics, it is important for students to be able to relate the directly measureable quantities like pressure loss in pipeline flow to entropy generation in real flows.

The main objectives of this article are: (1) to briefly review the background in fluid mechanics and thermodynamics related to pressure loss, mechanical energy dissipation, and entropy generation in real flows; (2) to demonstrate the link between mechanical energy dissipation and entropy generation in real flows; (3) to carry out experimental work to determine friction factor as a function of Reynolds number and entropy generation rate as a function of fluid velocity in flow of emulsion-type test fluids in different diameter pipelines; and (4) to explain mechanistically the cause of entropy generation in real flows.

#### **2. Background**

#### *2.1. Fluid Mechanics*

Consider flow of a fluid through a stationary control volume. The macroscopic balance of any entity (mass, momentum, energy) over the control volume can be expressed as:

$$Input-Output + \text{Generation} = \text{Accumulation} \tag{1}$$

The application of the entity balance equation, Equation (1), to mass gives the following integral equation [1,2]:

$$-\oint \rho \left(\vec{n} \cdot \vec{V}\right) dA = \iiint \frac{\partial \rho}{\partial t} d\mathcal{S} \tag{2}$$

where <sup>ρ</sup> is the fluid density, *<sup>n</sup>*<sup>ˆ</sup> is the unit outward normal to the control surface, <sup>→</sup> *V* is the fluid velocity vector, *A* is the control surface area, *t* is the time, ϑ is the volume of the control volume, the cyclic double integral is the surface integral over the entire control surface, and the cyclic triple integral is the volume integral over the entire control volume. It should be noted that the surface integral <sup>ρ</sup> - *n*ˆ· → *V dA* is the net outward flow of mass across the entire control surface. The volume integral ∂ρ <sup>∂</sup>*<sup>t</sup> d*ϑ is the rate of accumulation of mass within the entire control volume (assumed to be fixed and non-deforming). For a control volume with one inlet and one outlet, the macroscopic or integral mass balance equation, Equation (2), under steady state condition, reduces to:

$$
\rho\_1 A\_1 V\_1 = \rho\_2 A\_2 V\_2 = \dot{m} \tag{3}
$$

where *<sup>V</sup>* is the fluid velocity, *<sup>A</sup>* is the cross-section area of inlet or outlet opening, and . *m* is the mass flow rate. The subscript 1 indicates inlet variables and the subscript 2 indicates outlet variables.

The application of the entity balance equation, Equation (1), to linear momentum gives the following integral equation [1,2]:

$$-\oint \vec{\rho} \cdot \vec{\mathcal{V}} \left(\hat{n} \cdot \vec{\mathcal{V}}\right) dA + \sum \vec{F} = \iiint \oint \frac{\partial \left(\vec{\rho} \cdot \vec{\mathcal{V}}\right)}{\partial t} d\mathcal{S} \tag{4}$$

where <sup>→</sup> *F* is the force vector. Equation (4) is often referred to as *momentum theorem* [1]. The surface integral <sup>ρ</sup> → *V* - *n*ˆ· → *V dA* is the net outward flow of linear momentum across the entire control surface. The volume integral <sup>∂</sup> - ρ → *V* 

<sup>∂</sup>*<sup>t</sup> d*ϑ is the rate of accumulation of linear momentum within the entire control volume (assumed to be fixed and non-deforming). The momentum balance equation, Equation (4), is a vector equation. In Cartesian coordinates, the vector momentum balance equation can be written as three scalar equations:

$$-\iint \rho V\_{\mathbf{x}} \left(\hat{n} \cdot \overrightarrow{V}\right) dA + \sum F\_{\mathbf{x}} = \iiint \oint \frac{\partial (\rho V\_{\mathbf{x}})}{\partial t} dS \tag{5}$$

$$-\iint \rho \, V\_y \left(\hbar \vec{V}\right) dA + \sum F\_y = \iiint \frac{\partial \left(\rho V\_y\right)}{\partial t} d\mathcal{S} \tag{6}$$

$$-\iiint \rho V\_z \Big(\hat{n} \cdot \vec{V}\Big) dA + \sum F\_z = \iiint \oint \frac{\partial (\rho V\_z)}{\partial t} dS \tag{7}$$

where *Vx*, *Vy*, *Vz* and *Fx*, *Fy*, *Fz* are velocity and force components in *x*, *y*, *z* directions. For a control volume with one inlet and one outlet, the macroscopic or integral momentum balance equations, Equations (5)–(7), under steady state condition, reduce to:

$$V\_{\mathbf{x},1}(\rho\_1 A\_1 V\_1) - V\_{\mathbf{x},2}(\rho\_2 A\_2 V\_2) + \sum F\_{\mathbf{x}} = 0\tag{8}$$

$$V\_{y,1}(\rho\_1 A\_1 V\_1) - V\_{y,2}(\rho\_2 A\_2 V\_2) + \sum F\_y = 0\tag{9}$$

$$V\_{z,1}(\rho\_1 A\_1 V\_1) - V\_{z,2}(\rho\_2 A\_2 V\_2) + \sum F\_z = 0\tag{10}$$

Using mass balance, Equation (3), Equations (8)–(10) can be further simplified as:

$$\sum F\_{\mathbf{x}} = \dot{m}(V\_{\mathbf{x},2} - V\_{\mathbf{x},1}) \tag{11}$$

$$\sum F\_y = \dot{m} \left( V\_{y,2} - V\_{y,1} \right) \tag{12}$$

$$\sum F\_z = \dot{m}(V\_{z,2} - V\_{z,1})\tag{13}$$

Note that we are assuming velocity profiles are uniform at inlet and outlet of the control volume. The application of the entity balance equation, Equation (1), to mechanical energy gives the following integral equation assuming *frictionless* flow [1,2]:

$$-\iiint \rho \left(\not p \cdot \vec{V}\right) \left(\not p + KE + \frac{P}{\rho}\right) dA - \dot{W}\_{sh} = \iiint \not p \frac{\partial \left[\not p \left(\not p + KE\right)\right]}{\partial t} d\theta \tag{14}$$

where ϕ is the potential energy per unit mass, *KE* is the kinetic energy per unit mass, *P* is the pressure, and . *Wsh* is the rate of shaft work. The surface integral <sup>ρ</sup> - *n*ˆ· → *V* ϕ + *KE* + *<sup>P</sup>* ρ *dA* is the net outward flow of mechanical energy (including flow work) across the entire control surface. The volume integral <sup>∂</sup>[ρ(ϕ+*KE*)] <sup>∂</sup>*<sup>t</sup> d*ϑ is the rate of accumulation of mechanical energy within the entire control volume (assumed to be fixed and non-deforming). For a control volume with one inlet and one outlet, the macroscopic mechanical energy balance equation for frictionless flow, Equation (14), under steady state condition, reduces to:

$$\left(\varphi\_1 + KE\_1 + \frac{P\_1}{\rho\_1}\right) (\rho\_1 A\_1 V\_1) - \left(\varphi\_2 + KE\_2 + \frac{P\_2}{\rho\_2}\right) (\rho\_2 A\_2 V\_2) - \dot{W}\_{\text{sh}} = 0\tag{15}$$

Using mass balance, Equation (3), Equation (15) can be further simplified as:

$$
\dot{W}\_{sh} + \dot{m} \left[ \Delta q + \Delta (KE) + \Delta \left( \frac{P}{\rho} \right) \right] = 0 \tag{16}
$$

Note that mechanical energy is conserved only in frictionless flows. In real flows, however, mechanical energy is not conserved as mechanical energy dissipation occurs due to friction in fluids. The mechanical energy loss due to friction in real flows can be treated as negative generation (destruction) in entity balance equation, Equation (1). Thus, Equation (16) can be modified as [2]:

$$
\dot{\mathcal{W}}\_{sh} + \dot{F}\_l + \dot{m} \left[ \Delta q + \Delta (KE) + \Delta \left( \frac{P}{\rho} \right) \right] = 0 \tag{17}
$$

where . *Fl* is the rate of mechanical energy dissipation due to friction in fluid. Equation (17) assumes that the velocity profiles are uniform at inlet and outlet of the control volume. Furthermore any work due to viscous effects (shear stresses and viscous normal stresses) at the control surface is assumed to be negligible.

#### *2.2. Thermodynamics*

The *first law of thermodynamics* is simply the principle of conservation of energy. When applied to a control volume, it indicates that that the rate of accumulation of total energy inside the control volume is equal to the net rate of total energy addition to the control volume. Equation (1) is applicable to total energy as the entity with no generation. For a stationary control volume (fixed and non-deforming control volume), the first law of thermodynamics can be expressed as [1]:

$$-\iiint \rho \left(\vec{n} \cdot \vec{V}\right) \left(\varepsilon + \frac{P}{\rho}\right) dA + \dot{Q} - \dot{W}\_{sh} = \iiint \frac{\partial (\rho \varepsilon)}{\partial t} d\mathcal{S} \tag{18}$$

where *<sup>e</sup>* is the specific total energy (total energy per unit mass) of fluid and . *Q* is the rate of heat transfer. The specific total energy *e* includes internal energy, kinetic energy and potential energy, that is, *e* = *u* + *KE* + ϕ, where *u* is the specific internal energy. The work associated with shear stress and viscous portion of normal stress at the control surface is assumed to be zero [1]. The surface integral <sup>ρ</sup> - *n*ˆ· → *V e* + *<sup>P</sup>* ρ *dA* is the net outward flow of total energy (including flow work) across the entire control surface. The volume integral <sup>∂</sup>(ρ*e*) <sup>∂</sup>*<sup>t</sup> d*ϑ is the rate of accumulation of total energy within the entire control volume.

For a control volume with one inlet and one outlet, the macroscopic total energy balance (first law of thermodynamics), Equation (18), under steady state condition, reduces to:

$$(e\_1 + \frac{P\_1}{\rho\_1})(\rho\_1 A\_1 V\_1) - \left(e\_2 + \frac{P\_2}{\rho\_2}\right)(\rho\_2 A\_2 V\_2) + \dot{Q} - \dot{W}\_{sh} = 0\tag{19}$$

Using mass balance, Equation (3), Equation (19) can be further simplified as:

$$
\dot{m} \left[ \Delta \left( \varepsilon + \frac{P}{\rho} \right) \right] = \dot{Q} - \dot{W}\_{\text{sh}} \tag{20}
$$

As

$$\varepsilon + \frac{P}{\rho} = u + V^2/2 + gz + \frac{P}{\rho} = h + V^2/2 + gz \tag{21}$$

where *g* is acceleration due to gravity, *z* is elevation and *h* is specific enthalpy of fluid given as: *h* = *u* + *P*/ρ, Equation (20) could also be written as:

$$
\dot{m}\left[\Delta\text{h} + \Delta\left(V^2/2\right) + \text{g}\Delta z\right] = \dot{Q} - \dot{W}\_{sh} \tag{22}
$$

For reversible (frictionless) flows, it can be readily shown that

$$
\dot{Q} = \dot{m} \left[ \Delta h - \Delta \left( \frac{P}{\rho} \right) \right] \tag{23}
$$

Upon substitution of . *Q* from Equation (23) into Equation (22), the mechanical energy balance equation for frictionless flow, Equation (16), is recovered.

The *second law of thermodynamics* states that all irreversible (real) processes are accompanied by entropy generation in the universe [3]. For flow through a stationary control volume, the entropy generation rate ( . *SG*) in the universe can be expressed as:

$$\dot{S}\_{\text{G,universe}} = \dot{S}\_{\text{G,CV}} + \dot{S}\_{\text{G,Surr}} = \iiint \rho s \left(\hat{n} \cdot \overrightarrow{V}\right) dA + \iiint \frac{\partial(\rho s)}{\partial t} d\boldsymbol{\aleph} - \sum \frac{\dot{Q}\_i}{T\_i} \ge 0 \tag{24}$$

where *<sup>s</sup>* is the entropy per unit mass of fluid, . *Qi* is the rate of heat transfer to control volume from *ith* heat reservoir at an absolute temperature of *Ti*, the subscripts *CV* and *Surr* refer to control volume and surroundings, respectively. The equality in Equation (24) is valid for any reversible (frictionless) process and the inequality is valid for all irreversible processes. The surface integral <sup>ρ</sup>*<sup>s</sup>* - *n*ˆ· → *V dA* is the net outward flow of entropy across the entire control surface. The volume integral <sup>∂</sup>(ρ*s*) <sup>∂</sup>*<sup>t</sup> d*ϑ is the rate of accumulation of entropy within the entire control volume (assumed to be fixed and non-deforming).

For a control volume with one inlet and one outlet, Equation (24), under steady state condition, reduces to:

$$\dot{S}\_{\text{G,universe}} = (s\_2)(\rho\_2 A\_2 V\_2) - (s\_1)(\rho\_1 A\_1 V\_1) - \sum \frac{\dot{Q}\_i}{T\_i} \ge 0 \tag{25}$$

Using mass balance, Equation (3), Equation (25) can be further simplified as:

$$
\dot{S}\_{\text{G,universal}} = \dot{m}(\Delta \mathbf{s}) - \sum \frac{\dot{Q}\_i}{T\_i} \ge 0 \tag{26}
$$

#### *2.3. Steady Flow in a Pipe*

Consider steady flow of an incompressible fluid in a cylindrical pipe of uniform diameter (see Figure 1). From mass balance, .

$$
\dot{m} = \rho V\\\dot{A} = \text{constant} \tag{27}
$$

As ρ and *A* are constant, the velocity is constant,

$$V = \text{constant} \tag{28}$$

**Figure 1.** Stresses acting on fluid inside the control volume.

From momentum balance, Equation (11):

$$\sum F\_x = \dot{m}(V\_{x,2} - V\_{x,1}) = 0\tag{29}$$

Now consider force balance over a differential control volume of length *dx* as shown in Figure 1:

$$\sum F\_x = \begin{bmatrix} PA - (P + dP)A - \tau\_W(\pi D)dx = \end{bmatrix} \tag{30}$$

where τ*<sup>w</sup>* is the wall shear stress, *D* is the pipe internal diameter, and *A* = π*D*2/4. Equation (30) leads to:

$$-\left(dP/dx\right) = \frac{\pi D \tau\_w}{A} = \frac{4\tau\_w}{D} \tag{31}$$

Upon rearrangement, Equation (31) gives

$$
\tau\_{\rm av} = -\frac{(dP/dx)D}{4} = -\frac{(dP/dx)R}{2} \tag{32}
$$

where *R* is the pipe radius. Equation (32) can also be applied to any radial position as:

$$
\pi = -\frac{(dP/dx)r}{2} \tag{33}
$$

where *r* is any radial position in the pipe. From Equations (32) and (33), it follows that:

$$\frac{\pi}{\pi\_w} = \frac{r}{R} \tag{34}$$

Equation (34) describes the variation of shear stress with the radial position. The shear stress varies linearly with the radial position.

Using macroscopic mechanical energy balance, Equation (17), it can be readily shown that for steady incompressible flow in a horizontal pipe of uniform diameter:

$$
\dot{F}\_l/L = \left(\dot{m}/\rho\right)(-dP/dx)\tag{35}
$$

where *<sup>L</sup>* is the length of the pipe. The rate of mechanical energy loss . *Fl* can further be expressed in terms of a friction factor (*f*) defined as:

$$f = 2\tau\_w / \left(\rho \overline{V}^2\right) \tag{36}$$

where *V* is the average velocity in the pipe. Upon substitution of τ*<sup>w</sup>* from Equation (32) into Equation (36), we get

$$f = (-dP/dx)(D) / \left(2\rho \overline{V}^2\right) \tag{37}$$

From Equations (35) and (37), it follows that the mechanical energy loss per unit mass of fluid is:

$$
\dot{F}\_l/\dot{m} = 4f(L/D)\left(\overline{V}^2/2\right) \tag{38}
$$

From Equation (38), it follows that the mechanical energy loss per unit length per unit mass of fluid is: 

$$\dot{F}\_l' = \dot{F}\_l / (Lm) = (4f/D) \left(\overline{V}^2 / 2\right) \tag{39}$$

In order to calculate the mechanical energy loss in pipeline flows, the value of friction factor is required.

In *laminar* flow of Newtonian fluids (*Re* ≤ 2100), friction factor is related to Reynolds number through the following theoretical relationship [1,4]:

$$f = 16/\text{Re}$$

where the Reynolds number *Re* is defined as:

$$
\mathcal{R}\mathfrak{e} = \mathfrak{p}D\overline{V}/\mathfrak{\mu} \tag{41}
$$

In *turbulent* flow of Newtonian fluids, friction factor is a function of Reynolds number and relative roughness of pipe (/*D*). For hydraulically smooth pipes (/*D* → 0), the friction factor depends only on *Re* in turbulent regime. The following semi-empirical equation, often referred to as von Karman–Nikuradse equation, describes the *f vs*. *Re* turbulent behavior of Newtonian fluids in smooth pipes very well [1]:

$$1 \; \sqrt{f} = 4 \log\_{10} \left( \text{Re} \, \sqrt{f} \right) - 0.40 \tag{42}$$

The von Karman–Nikuradse equation is not explicit in friction factor. A number of explicit *f vs*. *Re* relations are available in the literature. One of the popular ones is the Blasius friction factor equation for turbulent flow of Newtonian fluids in smooth pipes [4]:

$$f = 0.079 / \text{R}\omega^{0.25} \tag{43}$$

Equation (43) is accurate over a Reynolds number range of 3000 ≤ *Re* ≤ 100,000. For turbulent flow in rough pipes, the following Colebrook equation [5,6] is widely accepted:

$$1 \,\,\sqrt{f} = -4\log\_{10}\left(\frac{\varepsilon/D}{3.7} + \frac{1.26}{\text{Re}\,\sqrt{f}}\right) \tag{44}$$

This Colebrook equation is implicit in friction factor. A number of explicit *f vs*. *Re* relations are available in the literature for turbulent flow of Newtonian fluids in rough pipes [7–9]. An explicit equation which is very accurate for turbulent flow of Newtonian fluids in rough pipes is as follows:

$$1/\sqrt{f} = -4\log\_{10}\left[\frac{\epsilon/D}{3.7} - \frac{5.02}{Re}\log\_{10}\left(\frac{\epsilon/D}{3.7} - \frac{5.02A}{Re}\right)\right] \tag{45}$$

$$A = \log\_{10} \left( \frac{\epsilon / D}{3 .7} + \frac{13}{Re} \right) \tag{46}$$

This equation was originally proposed by Zigrang and Sylvester [10].

#### *2.4. Entropy Generation in Steady Flow in a Pipe*

In flow of real fluids, the dissipation of mechanical energy, and hence loss of pressure, is simultaneously reflected in the generation of entropy [11,12]. Consequently, the pipeline flow experiments performed by undergraduate students in the fluid mechanics laboratory can also be used as a tool to teach the second law of thermodynamics which states that all real processes are accompanied by generation of entropy in the universe.

For steady flow in a pipe with no heat transfer, Equation (26) reduces to:

$$
\dot{S}\_{\text{G,universe}} = \dot{S}\_{\text{G,CV}} = \dot{m}(\Delta \text{s}) > 0 \tag{47}
$$

There is no entropy generation in the surroundings. All the entropy is generated within the fluid inside the pipe and the rate of entropy generation is the net rate of increase in entropy of the flowing stream.

We can now relate entropy change of the fluid stream to other variables. For pure substances, the relationship between entropy and other state variables is given as [3]:

$$Tds = dh - \left(dP/\rho\right)\tag{48}$$

where *T* is the absolute temperature. From the first law of thermodynamics, Equation (22), the enthalpy change is zero in the absence of heat transfer and shaft work for steady flow in a horizontal pipe of uniform diameter. Consequently, Equation (48) reduces to:

$$T\text{ds} = -(dP/\rho)\tag{49}$$

Assuming incompressible flow and constant temperature, Equation (49) upon integration gives:

$$
\Delta \mathbf{s} = -\frac{\Delta P}{\rho T} \tag{50}
$$

Strictly speaking, the temperature is expected to rise somewhat in adiabatic flow due to frictional heating. However, the temperature rise is usually very small in pipeline flow experiments conducted in the undergraduate fluid mechanics laboratory. From Equations (47) and (50), it follows that:

$$
\dot{S}\_G = \frac{\dot{m}}{\rho} \left( -\frac{\Delta P}{T} \right) > 0 \tag{51}
$$

The subscript "*CV*" has been removed from . *SG*,*CV* for the sake of simplicity. We can also express the rate of entropy generation in a pipe on a unit length basis as:

$$\dot{S}\_G^{'} = \dot{m} \left( \frac{F\_l^{'}}{T} \right) = \frac{\dot{m}}{\rho T} \left( -\frac{dP}{d\mathbf{x}} \right) \tag{52}$$

where . *S <sup>G</sup>* is the rate of entropy generation per unit length of the pipe. From Equations (37) and (52), it can be readily shown that: .

$$
\dot{S}\_G' = \left(\frac{\pi}{2T}\right) \left(\rho D \overline{V}^3\right) \tag{53}
$$

In laminar flow, the friction factor *f* is given by Equation (40). Consequently, Equation (53) yields:

$$
\dot{S}\_G' = \left(\frac{8\pi}{T}\right)\mu \overline{V}^2\tag{54}
$$

Thus entropy generation rate per unit length of pipe in steady laminar flow of a Newtonian fluid is directly proportional to fluid viscosity and square of average velocity in the pipe.

In turbulent flow of a Newtonian fluid in hydraulically smooth pipe, Equations (43) and (53) give the flowing expression for entropy generation rate per unit length of pipe:

$$\dot{S}\_G^{'} = \left(\frac{0.079\pi}{2T}\right) \mu^{0.25} \rho^{0.75} D^{0.75} \overline{V}^{2.75} \tag{55}$$

In turbulent flow, the entropy generation rate per unit length of pipe also depends on pipe diameter and fluid density, in addition to viscosity and fluid velocity dependence. Although the viscosity dependence of . *S <sup>G</sup>* in turbulent flow is less severe in comparison with laminar flow, the velocity dependence is stronger in turbulent flow.

Figure <sup>2</sup> shows the plots of . *S <sup>G</sup> vs*. *V* on a log-log scale for laminar and turbulent flows generated from Equations (54) and (55), respectively. The fluid properties used in the equations are: μ = 1 mPa·s, <sup>ρ</sup> = 1000 kg/m3. The temperature used is 298.15 K. A single line of slope 2 is obtained for laminar regime regardless of the pipe diameter. The entropy generation rate per unit length of the pipe increases linearly with the increase in average velocity in the pipe. A family of parallel lines of slope 2.75 is obtained for the turbulent regime. The line shifts upward towards higher entropy generation rate with the increase in the pipe diameter.

**Figure 2.** Entropy generation rate per unit length of pipe as a function of average velocity in the pipe.

#### **3. Experimental Work**

#### *3.1. Apparatus*

A flow rig consisting of five different diameter pipeline test sections (stainless steel seamless tubes, hydraulically smooth) was designed and constructed. The pipelines were installed horizontally. Table 1 gives the dimensions of the test sections. The test fluid was circulated through the pipeline test sections, one at a time, using a centrifugal pump. An electromagnetic flow meter was used to measure the flow rate of a fluid circulated through the test section. The pressure drop in a pipeline test section was measured using pressure transducers covering a broad range of pressure drops. The pressure drop as a function of flow rate was recorded by a computer data acquisition system. The experiments were carried out at a constant temperature of 25 ◦C.


**Table 1.** Various dimensions of the pipeline flow test sections.

#### *3.2. Test Fluids*

The test fluids used were surfactant-stabilized oil-in-water (O/W) emulsions. The viscosity of the test fluid was increased by increasing the oil concentration of the emulsion. Table 2 summarizes the viscosity and density data of the test fluids at 25 ◦C.


**Table 2.** Viscosity and density of test fluids.

#### **4. Results and Discussion**

#### *4.1. Fluid Mechanics Experiments*

The experimental data obtained from pipeline flow experiments consist of pressure drop (over a known length of pipe) versus flow rate for different diameter pipes. From pressure drop versus flow rate data, friction factor is calculated from Equation (37), re-written as:

$$f = (-\Delta P/L)(D) / \left(2\rho \overline{V}^2\right) \tag{56}$$

where the average velocity *V* is obtained from:

$$
\overline{V} = 4\dot{\Phi}/\left(\pi D^2\right) \tag{57}
$$

Here . ϑ is the volumetric flow rate of fluid. The Reynolds is calculated from the defining relation, Equation (41), re-written as: .

$$
\mathcal{Re} = 4\rho\theta / (\pi D \mu) \tag{58}
$$

The *f vs*. *Re* data are plotted on a log-log scale and compared with the predictions of available *f vs*. *Re* theoretical and empirical relations (Equation (40) for laminar flow and Equation (43) for turbulent flow in hydraulically smooth pipes).

Figures 3–8 show the *f vs*. *Re* data for different test fluids obtained from different diameter pipelines. The experimental data follow the existing *f vs*. *Re* relations reasonably well in both laminar and turbulent regimes.

**Figure 3.** Friction factor versus Reynolds number data for aqueous surfactant solution.

**Figure 4.** Friction factor versus Reynolds number data for 16.53% by volume oil-in-water emulsion.

**Figure 5.** Friction factor versus Reynolds number data for 30.4% by volume oil-in-water emulsion.

**Figure 6.** Friction factor versus Reynolds number data for 44.41% by volume oil-in-water emulsion.

**Figure 7.** Friction factor versus Reynolds number data for 49.65% by volume oil-in-water emulsion.

**Figure 8.** Friction factor versus Reynolds number data for 55.14% by volume oil-in-water emulsion.

#### *4.2. Entropy Generation Results*

The entropy generation rate in pipeline flow is calculated from pressure drop versus flow rate measurements using Equation (52), re-written as:

$$
\dot{S}\_G' = \frac{\dot{\Phi}}{T} \left( -\frac{\Delta P}{L} \right) \tag{59}
$$

where *T* = 25 + 273.15 = 298.15*K*.

The . *S <sup>G</sup> vs*. *V* data are plotted on a log-log scale and compared with the predictions of equations developed in Section 2.4, that is, Equation (54) for laminar flows and Equation (55) for turbulent flows in smooth pipes.

Figures 9–14 show the . *S <sup>G</sup> vs*. *V* data for different test fluids obtained from different diameter pipelines. As expected from Equation (54), the experimental data corresponding to laminar flow is independent of the pipe diameter. Also, the slope of the laminar flow line is 2 as predicted by Equation (54). Equation (54) describes the laminar flow . *S <sup>G</sup> vs*. *V* data adequately for all the test fluids investigated. According to Equation (55), the turbulent flow . *S <sup>G</sup> vs*. *V* data for a given diameter pipeline

should follow a straight line on a log-log scale with a slope of 2.75. As expected, the experimental data in turbulent regime follow a straight line of slope 2.75. With the increase in pipe diameter, the entropy generation rate increases. The . *S <sup>G</sup> vs*. *V* line shifts upward with the increase in pipe diameter but the slope remains the same, that is, 2.75. Thus Equation (55) describes entropy generation in turbulent flows reasonably well.

**Figure 9.** Entropy generation rate per unit length of pipe as a function of average velocity in flow of aqueous surfactant solution.

**Figure 10.** Entropy generation rate per unit length of pipe as a function of average velocity in flow of 16.53% by volume oil-in-water emulsion.

**Figure 11.** Entropy generation rate per unit length of pipe as a function of average velocity in flow of 30.4% by volume oil-in-water emulsion.

**Figure 12.** Entropy generation rate per unit length of pipe as a function of average velocity in flow of 44.41% by volume oil-in-water emulsion.

**Figure 13.** Entropy generation rate per unit length of pipe as a function of average velocity in flow of 49.65% by volume oil-in-water emulsion.

**Figure 14.** Entropy generation rate per unit length of pipe as a function of average velocity in flow of 55.14% by volume oil-in-water emulsion.

#### *4.3. Mechanism of Entropy Generation in Real Flows*

In real flows, entropy generation is a volumetric phenomenon caused by friction (non-zero viscosity) in fluids. Due to non-zero viscosity of real fluids, velocity gradients and viscous stresses are set up when fluid is forced to flow through a pipe. The presence of viscous stresses and velocity gradients in the fluid cause mechanical energy dissipation into frictional heating effect (internal energy). The degradation of highly ordered mechanical energy into disorderly internal energy results in entropy generation.

The entropy generation rate per unit volume of fluid in real flows ( . *S <sup>G</sup>* ) can be expressed as [13]:

$$\vec{S}\_{G}^{\prime\prime\prime} = \frac{1}{T} \begin{pmatrix} \underline{\tau} & : \vec{\nabla} \vec{V} \end{pmatrix} \tag{60}$$

where τ <sup>=</sup> is the viscous stress tensor and <sup>∇</sup> → *V* is the velocity gradient tensor. It should be noted that we are assuming that there are no temperature gradients in the fluid. If temperature is not uniform and temperature gradients are present, then entropy generation can occur due to another mechanism, namely, irreversible transfer of heat.

In general, the local entropy generation rate per unit volume of fluid ( . *S <sup>G</sup>* ) will vary with the position coordinates. As an example, we illustrate the application of Equation (60) to steady laminar flow in cylindrical pipes. For steady laminar flow of a fluid in a uniform diameter pipe, the term τ <sup>=</sup> : <sup>∇</sup> → *V* simplifies to:

$$
\tau\_{\equiv} : \vec{\nabla V} = \tau\_{r\chi} \left( \frac{\partial V\_{\chi}}{\partial r} \right) \tag{61}
$$

From Equations (60) and (61), it follows that:

$$\vec{S}\_{G}^{\prime\prime\prime} = \frac{1}{T} \begin{pmatrix} \tau & : \vec{V} \\ \underline{\phantom{\prime}} & : \vec{V} \end{pmatrix} = \frac{1}{T} \tau\_{\text{rx}} \begin{pmatrix} \partial V\_{\text{x}} \\ \partial \boldsymbol{r} \end{pmatrix} \tag{62}$$

From Newton's law of viscosity,

$$
\pi\_{rx} = \mu \left(\frac{\partial V\_x}{\partial r}\right) \tag{63}
$$

Consequently,

$$\dot{S}\_{G}^{\prime\prime\prime} = \frac{1}{T} \pi\_{r\mathbf{x}} \left(\frac{\partial V\_{\mathbf{x}}}{\partial r}\right) = \frac{\mu}{T} \left(\frac{\partial V\_{\mathbf{x}}}{\partial r}\right)^{2} \tag{64}$$

The velocity distribution in steady laminar flow of a Newtonian fluid is given as:

$$V\_X = 2\overline{V} \left[ 1 - \left(\frac{r}{R}\right)^2 \right] \tag{65}$$

Thus,

$$\frac{\partial V\_x}{\partial r} = -\frac{4r}{R^2}\overline{V} \tag{66}$$

From Equations (64) and (66), we get:

$$\dot{S}\_{G}^{\prime\prime\prime} = \frac{\mu}{T} \left(\frac{\partial V\_{x}}{\partial r}\right)^{2} = \left(\frac{16\mu}{T}\right) \left(\frac{r^{2}}{R^{4}}\right) \overline{V}^{2} \tag{67}$$

This equation describes the local rate of entropy generation per unit volume in steady laminar flow of a Newtonian fluid in cylindrical pipe of uniform diameter. Figure <sup>15</sup> shows the plot of . *S <sup>G</sup>* as a function of radial position. The plot is generated using Equation (67). . *S <sup>G</sup>* is zero at the centre of the pipe as velocity gradiant and hence local dissipation of mechanical energy into frictional heating is zero at the centre of the pipe. As we go towards the pipe wall, the local entropy generation increases due to an increase in the velocity gradient and hence frictional heating of the fluid.

**Figure 15.** Local entropy generation rate per unit volume ( . *S <sup>G</sup>* ) as a function of radial position in laminar flow of a Newtonian fluid (*R* = 10 mm, *T* = 298.15 K, μ = 10 mPa·s, *V* = 0.1 m/s).

Finally it can be readily shown that the global rate of entropy generation (Equation (54)) follows from the integration of the local rate of entropy generation. The global entropy generation rate per unit length . *S <sup>G</sup>* can be expressed as:

$$\dot{S}\_G^{'} = \frac{1}{T} \int\_0^R \left(\dot{S}\_G^{\prime\prime\prime}\right) (2\pi r) dr\tag{68}$$

From Equations (67) and (68), we get:

$$\dot{S}\_G' = \frac{32\pi\mu}{TR^4} \overline{\mathcal{V}}^2 \int\_0^R r^3 dr = \left(\frac{8\pi}{T}\right) \mu \overline{\mathcal{V}}^2\tag{69}$$

This is the same result that was obtained earlier in Equation (54).

#### **5. Conclusions**

In conclusion, a novel approach is described to teach the second law of thermodynamics with the help of an undergraduate fluid mechanics laboratory involving pipeline flow experiments. The relevant background in fluid mechanics and thermodynamics is reviewed briefly. The link between entropy generation and pressure loss in pipeline flow experiments is demonstrated both theoretically and experimentally. Experimental work involving flow of emulsion-type test fluids in different diameter pipes is carried out to determine friction factor versus Reynolds number behavior and entropy generation rates in pipeline flows. Entropy generation in pipeline flows is explained mechanistically considering local entropy generation in the presence of viscous stresses and velocity gradients in flow of real fluids.

**Funding:** This research received no external funding. **Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **What Can Students Learn While Solving Colebrook's Flow Friction Equation?**

**Dejan Brki´c 1,2,\*,**† **and Pavel Praks 1,\*,**†


Received: 12 May 2019; Accepted: 25 June 2019; Published: 27 June 2019

**Abstract:** Even a relatively simple equation such as Colebrook's offers a lot of possibilities to students to increase their computational skills. The Colebrook's equation is implicit in the flow friction factor and, therefore, it needs to be solved iteratively or using explicit approximations, which need to be developed using different approaches. Various procedures can be used for iterative methods, such as single the fixed-point iterative method, Newton–Raphson, and other types of multi-point iterative methods, iterative methods in a combination with Padé polynomials, special functions such as Lambert W, artificial intelligence such as neural networks, etc. In addition, to develop explicit approximations or to improve their accuracy, regression analysis, genetic algorithms, and curve fitting techniques can be used too. In this learning numerical exercise, a few numerical examples will be shown along with the explanation of the estimated pedagogical impact for university students. Students can see what the difference is between the classical vs. floating-point algebra used in computers.

**Keywords:** Colebrook equation; Lambert W function; Padé polynomials; iterative methods; explicit approximations; learning; teaching strategies; floating-point computations

#### **1. Introduction**

The Colebrook equation for flow friction is an empirical formula developed by C. F. Colebrook in 1939 [1] and it is based on his joint experiment, which was performed with Prof. C. M. White in 1937 [2]. Colebrook was at this time a PhD student of Prof. White at Imperial College in London, UK [3]. The experiment was performed with a set of pipes of which the inner side was covered with sand of different grain sizes, while one was left completely without sand to be smooth. The sand was fixed on the pipe's inner surface using a sort of glue as adhesive material, while air was used as the fluid in this experiment.

The Colebrook equation follows logarithmic law and, at first sight, it seems to be very simple. However, it is given implicitly, with respect to the unknown flow friction factor *f* in a way that it can be expressed explicitly only in terms of the Lambert W function, with further difficulties in its evaluation [4–6] (it can be rearranged for gas flow [7]). Otherwise, it needs to be solved iteratively [8–10] or explicitly, using approximate formulas specially developed for such a purpose [10–17]. The unknown flow friction factor *f* in the Colebrook's equation depends on the known Reynolds number *Re* and the relative roughness of the inner pipe surface <sup>ε</sup> *<sup>D</sup>* , all three of which are dimensionless, shown in Equation (1):

$$\frac{1}{\sqrt{f}} = -2 \cdot \log\_{10} \left( \frac{2.51}{\text{Re}} \cdot \frac{1}{\sqrt{f}} + \frac{\frac{\varepsilon}{D}}{3.71} \right) \tag{1}$$

The relative roughness <sup>ε</sup> *<sup>D</sup>* can take values from 0 to 0.05, while the Reynolds number *Re* is from about 2500 to 10<sup>8</sup> [18]. The Colebrook equation is valid for the turbulent condition of flow [19].

In this learning exercise, students can have numerous numerical tasks that can improve their computational skills. The task can start with a simple trial/error approach that can, in a certain way, put the Colebrook equation in balance, while, furthermore, a simple fixed-point iteration method can be introduced in the curriculum [20]. Then, more complex approaches that require derivatives of the Colebrook equation can lead students to the Newton–Raphson iterative methods [21–23] or to the more complex multi-point iterative methods [8,9]. A first simple explicit approximate formula with inner iterative cycles can be introduced using such approaches [24].

Transcendental functions, such as logarithmic and exponential functions, on computers require the execution of additional floating point operations [25–27]. On the contrary, basic mathematical operations (addition, subtraction, multiplication, and division) are very fast on computers, because they are executed directly. Students can learn how to avoid using these time-consuming functions [27–30]. This simplification process can have a large effect on the students' future jobs, regarding the design of complex pipe networks, which involve extensive computations [31–36]. For example, in Section 2.2.1, students can reuse the value of the logarithm, which is computed only in the first iteration, to all subsequent iterations [30]. The less computationally demanded Padé polynomials are used instead of the logarithmic function in all subsequent iterations [30]. The students in that way need to know how to generate the appropriate polynomial expression (using appropriate software tools).

Regarding special functions, the Lambert W function can be introduced in the curriculum [5,6], as the Lambert W function transforms the Colebrook equation in the explicit form [4]. The widely used explicit form has one fast-growing term that makes calculations impossible in computers, due to an overflow error [37,38]. Students should figure out why such a numerical error appears in computers and then understand how such situations can be avoided. For this reason, the next step involves the series expansion of the Lambert W function, and its cognate, the Wright ω function [39]. In this step, students will learn how to make accurate, but at the same time, efficient approximations of the Colebrook equation using the Wright ω function [27–29,40,41]. Last, but not least, students can use raw numerical data that can efficiently simulate the Colebrook equation and, in that way, learn principles of using artificial intelligence in everyday engineering practice [12,42–45].

Finally, keeping in mind that the Moody diagram represents a graphical interpretation of the Colebrook equation, students should make a sketch of the Moody diagram [18,46], which will add an additional value both from the engineering and mathematical points of view (Figure 1 of this article can be used for this purpose).

#### **2. Tasks for Students**

#### *2.1. Iterative Methods*

For the purpose of this case study, two numerical examples are given: (1) *Re* <sup>=</sup> 2.3 <sup>×</sup> <sup>105</sup> and <sup>ε</sup> *<sup>D</sup>* <sup>=</sup> <sup>10</sup><sup>−</sup>4; (2) *Re* <sup>=</sup> 4.6 <sup>×</sup> 107 and <sup>ε</sup> *<sup>D</sup>* = 0.037. Corresponding solutions for the flow friction factor are (1) *f* = 0.016050961, and (2) *f* = 0.062427396, respectively. Positions of these two examples with respect to the hydraulic regimes (laminar, transition to turbulent, and turbulent, which includes smooth, smooth to rough, and fully turbulent rough regime) are shown in Figure 1. Cases that can capture the nonlinearity of the Colebrook equation are: (a) Smooth pipe at very low, moderate, and very high *Re*, (b) pipe with moderate roughness at very low, moderate, and very high *Re*, and (c) pipe with high roughness ( <sup>ε</sup> *<sup>D</sup>* ∼ 0.05) and again, low, moderate, and high *Re*. Example 1 is in the zone of moderate values of the Reynolds number *Re*, where the proposed value of relative roughness, which is small ( ε *<sup>D</sup>* <sup>∼</sup> <sup>10</sup><sup>−</sup>4), is large enough to provide a transition from smooth towards the turbulent regime. On the other hand, Example 2 with higher values of *Re* and <sup>ε</sup> *<sup>D</sup>* is in a fully turbulent zone. Using Figure 1, other different examples can be provided for these teaching exercises.

*Fluids* **2019**, *4*, 114

Example 1 is set to be exactly in the zone where the smooth nature of the viscose zone of the inner pipe surface starts to be too thin and when roughness starts to have influence [47,48].

**Figure 1.** Position of the two presented examples with respect to the hydraulic regimes (laminar, transition to turbulent, and turbulent, which includes smooth, smooth to rough, and fully turbulent rough regime).

From Figure 1, students should see that our numerical Example 1 is set to be in the turbulent zone. More precisely, Example 1 is set at the end of the smooth turbulent regime, where the transition to fully turbulent begins. On the contrary, Example 2 is located in the fully developed rough turbulent zone. Students should understand the physical interpretation of different hydraulic regimes.

For this task, 45 min should be assigned.

#### 2.1.1. Trial/Error Method

To do this task, students should guess somehow the initial value of the friction factor *f*, then insert it into the Colebrook equation, and finally, calculate a new value of *f* until an error criterion is fulfilled. Using this iterative procedure, students should show how many iterations will probably be required, as the number of iterations depends on the initial guess and whether such a rudimentary iterative procedure will always lead to the final balanced solution or not. Here will be introduced a concept of convergence, divergence, and oscillations in iterative procedures. Students should understand how significant it is to find the adequate value for the starting point for an iterative procedure [9]. For our Example 1, which uses *Re* = 2.3 <sup>×</sup> 105 and <sup>ε</sup> *<sup>D</sup>* = <sup>10</sup>−4, the starting point needs to be chosen between −2.469 and 1574003 when Microsoft (MS) Excel is used and students should discuss why [10,20]. Subsequently, students should analyze the starting point for Example 2, in which *Re* <sup>=</sup> 4.6 <sup>×</sup> <sup>10</sup><sup>7</sup> and <sup>ε</sup> *<sup>D</sup>* = 0.037. Students should explain consequences of the inappropriately chosen initial starting point. Finally, students should understand that the starting point for −2.5, for example, will actually be

√1 *f*0 <sup>=</sup> <sup>ε</sup> *D* 3.71 <sup>−</sup> 2.5 <sup>=</sup> <sup>10</sup>−<sup>4</sup> 3.71 <sup>−</sup> 2.5 <sup>=</sup> <sup>−</sup>2.499973046 for our Example 1, in which <sup>ε</sup> *<sup>D</sup>* = <sup>10</sup>−<sup>4</sup> (the negative starting point does not have physical meaning and it is chosen purely for testing methods and to introduce the concept of "thinking outside of the box", which could contribute to development of creative thinking in general.

For this task, 45 min should be assigned.

#### 2.1.2. Fixed-Point Iterative Procedure

The fixed-point iterative method is suitable for spreadsheet solvers, such as MS Excel [10,20]. To prepare this task, students should rearrange the Colebrook equation in the suitable form *F*(*x*) = 0; where *x* = √<sup>1</sup> *f* , i.e., *<sup>x</sup>* <sup>+</sup> <sup>2</sup>· log10- 2.51 *Re* ·*<sup>x</sup>* <sup>+</sup> <sup>ε</sup> *D* 3.71 = 0. The iterative process starts with a selected starting point *x*<sup>0</sup> and continues by *xi*+<sup>1</sup> = *xi* − *F*(*xi*) until a stopping criteria is fulfilled, for example, when the results from two subsequent iterations are almost equal (*xi*+<sup>1</sup> − *xi*+<sup>1</sup> ∼ 0). From this approach, students will learn that the trial/error method can be used in a more systematic way. This task can improve the ability of students to use spreadsheet solvers, because students should set the solver to allow iterative computations. Moreover, the maximal number of iterations and the stopping criteria represented by the minimal value between two successive iterations should be set in the solver. For example, the MS Excel spreadsheet solver allows a maximal number of iterations of 32,767, and students should figure out why this constraint is set. The answer is that 215 = 32,768 as this represents the maximum it can be set for 16-bit registers in binary logic because one bit is reserved for the sign (±). In that way, students will learn that for human perception, more common is the decimal system, which is used, for example, for banknotes, \$1, \$10, \$100, while, for computers, the binary system is more suitable, which commonly uses 2, 4, 8, 16, 32, 64, etc.

In our case, (1) *Re* = 2.3 <sup>×</sup> 105 and <sup>ε</sup> *<sup>D</sup>* <sup>=</sup> <sup>10</sup>−4, for <sup>√</sup><sup>1</sup> *f*0 <sup>=</sup> <sup>ε</sup> *D* 3.71 − 2.4 ≈ −2.4 and for <sup>=</sup> <sup>ε</sup>

√1 *f*0 *D* 3.71 + 1,574,003 ≈ 1,574,003, the number of required iterations for nine decimal accuracy will be 11 in both cases, i.e., the accurate solution will be √<sup>1</sup> *f*11 = 7.8931. Students should make further tests with the same example, but with different starting points and, in addition, perform a similar test for Example 2.

Using the spreadsheet solver, students should first solve the Colebrook equation iteratively, using only one cell of the Colebrook equation for smooth conditions, when <sup>ε</sup> *<sup>D</sup>* is zero. Students should realize how important the role of the initial starting point is, because, in this case, software will declare a division by zero error. In other conditions, when <sup>ε</sup> *<sup>D</sup>* - 0, the initial starting point will be automatically set as *<sup>f</sup>*<sup>0</sup> <sup>=</sup> log10 ε *D* 3.71 .

Here, students may try to make a sketch of the function, i.e., to make a two-dimensional (2D) diagram of the Colebrook equation. Here, two possibilities should be analyzed: √<sup>1</sup> *<sup>f</sup>* <sup>=</sup> *<sup>F</sup> Re*· *f*, <sup>ε</sup> *D* 

as suggested by Rose [46], or *f* = *F Re*, <sup>ε</sup> *D* as suggested by Moody [18] (approach from Figure 1 of this article can be used). In addition, students can analyze whether a simple transformation *x* = √<sup>1</sup> *f*

accelerates solving of the Colebrook equation.

For this task, 45 min should be assigned.

2.1.3. Newton–Raphson Iterative Procedure

In this part, students will learn how to find the first derivative of the function with respect to the variable flow friction factor *f*, where *x* = √<sup>1</sup> *f* and *<sup>F</sup>*(*x*) <sup>=</sup> *<sup>x</sup>* <sup>+</sup> 2 log10- 2.51 *Re* ·*<sup>x</sup>* <sup>+</sup> <sup>ε</sup> *D* 3.71 , as in Equation (2):

$$\frac{\partial (F(\mathbf{x}))}{\partial \mathbf{x}} = F'(\mathbf{x}) = F'\left(\mathbf{x} + 2\log\_{10}\left(\frac{2.51}{Re}\cdot\mathbf{x} + \frac{\frac{\delta}{D}}{3.71}\right)\right) = -\frac{1}{2}\left(\frac{1}{\sqrt{|f|}}\right)^3 \cdot \left(1 + 2\cdot\frac{2.51}{Re \cdot \ln(10) \cdot \left(\frac{2.51}{3E}\cdot\frac{1}{\sqrt{|f|}} + \frac{\frac{\delta}{D}}{\sqrt{|f|}}\right)}\right) \tag{2}$$

The procedure for finding the first derivative can be done manually, but also for Equation (2), it can be generated using appropriate software, e.g., Matlab by the "diff" command.

The first derivative is used in the Newton–Raphson procedure as *fi*+<sup>1</sup> <sup>=</sup> *fi* <sup>−</sup> *<sup>F</sup>*(*fi*) *<sup>F</sup>*(*fi*) and students may test whether the Newton–Raphson procedure is faster or not, compared with the simple fixed-point method of the Colebrook equation, where *F*(*fi*+1) = √<sup>1</sup> *<sup>f</sup>* <sup>+</sup> 2 log10 2.51 *Re* · <sup>√</sup><sup>1</sup> *fi* <sup>+</sup> <sup>ε</sup> *D* 3.71 .

Students should realize that the previous example with the fixed point method is set for solving *x* = √<sup>1</sup> *f* , while here, the Newton–Raphson method is set for solving *F*(*x*) = 0. For Example 1, in which *Re* = 2.3 <sup>×</sup> <sup>10</sup><sup>5</sup> and <sup>ε</sup> *<sup>D</sup>* <sup>=</sup> <sup>10</sup><sup>−</sup>4, we can use the already analyzed starting point <sup>√</sup><sup>1</sup> *f*0 <sup>=</sup> <sup>ε</sup> *D* 3.71 + 1,574,003 ≈ 1,574,003 <sup>→</sup> *<sup>f</sup>*<sup>0</sup> = 4.04 <sup>×</sup> 10<sup>−</sup>13, which will give accurate results. Unfortunatelly, this iteration process requires 27 iterations, *f*<sup>27</sup> = 0.016050961, whereas √<sup>1</sup> *f*0 <sup>=</sup> <sup>ε</sup> *D* 3.71 − 2.4 ≈ −2.4 → *f*<sup>0</sup> = 0.173611 will report a division by zero error in the fourteenth iteration and the calculation will stop without obtaining the precise solution. Here, students should understand that even a more powerful iterative procedure can fail when the problem is not set in an appropriate way and when the role of a good initial starting point is not fully understood.

The Newton–Raphson method, in which *F* (*fi*) = 1 is assumed, is equivalent to the fixed-point method in the form *fi*+<sup>1</sup> = *fi* −*F*(*fi*). However, this formulation shows very bad convergence properties. Students should understand why the iterative process *xi*+<sup>1</sup> = √<sup>1</sup> *fi*+<sup>1</sup> = *xi* − *F*(*xi*) gives the fast solution, while the iterative process *fi*+<sup>1</sup> = *fi* − *F*(*fi*) is slow. Students should again analyze coordinates in diagrams from Rose and Moody, and discuss which formulation is better [18]. 

Here, it should be noted that the function *<sup>F</sup>*(*x*) <sup>=</sup> *<sup>x</sup>* <sup>+</sup> <sup>2</sup> log10-2.51 *Re* ·*<sup>x</sup>* <sup>+</sup> <sup>ε</sup> *D* 3.71 in form *F* √1 *f* = 2 log10 2.51 *Re* · <sup>√</sup><sup>1</sup> *<sup>f</sup>* <sup>+</sup> <sup>ε</sup> *D* 3.71 −<sup>2</sup> + √<sup>1</sup> *f* shows very slow convergence when it is used in the iterative process *fi*+<sup>1</sup> <sup>=</sup> *fi* <sup>−</sup> *<sup>F</sup>*(*fi*)

*<sup>F</sup>*(*fi*). Students should discover why. Students should investigate how the initial starting point affects the required number of iterations.

For this task, 90 min should be assigned (including repetition and clarifications from the previous tasks).

#### 2.1.4. Multi-Point Iterative Procedures

The shown Newton–Raphson method belongs to the two-point iterative methods and, in this task, students should find more powerful iterative methods. For example, the proposed methods for this exercise are three-point iterative methods, such as Halley or Schröder or even more powerful multi-point methods [8,9]. Students should figure out whether the number of required iterations to find the accurate solution decreases with the increased complexity of the chosen method. If the number of required iterations does not decrease, students should analyze why. In addition, it needs to be found how important it is to find the adequate initial starting point. In this task, students should find the second derivative of the Colebrook equation, which will increase their computing skills. In this exercise, students should develop the ability to make a balance between complexity, speed, and required accuracy. Students obtain an experience for choosing the appropriate iterative method.

Here, we will not repeat the Newton–Raphson calculations in the form of Equation (3), which gives fast accurate results, but we will show the Halley method: Equation (4) for which the first derivative *F* (*x*) and the second *F*(*x*), with respect to variable *x* = √<sup>1</sup> *fi* , are given by Equations (5) and (6), respectively.

$$\mathbf{x}\_{i+1} = \mathbf{x}\_i - \frac{F(\mathbf{x}\_i)}{F'(\mathbf{x}\_i)},\tag{3}$$

$$\mathbf{x}\_{i+1} = \mathbf{x}\_i - \frac{2 \cdot F(\mathbf{x}\_i) \cdot F'(\mathbf{x}\_i)}{2 \cdot F'(\mathbf{x}\_i)^2 - F(\mathbf{x}\_i) \cdot F''(\mathbf{x}\_i)} \,' \tag{4}$$

$$F'(\mathbf{x}\_i) = \frac{9287 \cdot \ln(10) \cdot \mathbf{x}\_i + 1000 \cdot \frac{\varepsilon}{D} \cdot \mathrm{Re} \cdot \ln(10) + 18574}{\ln(10) \cdot \left(9287 \cdot \mathbf{x}\_i + 1000 \cdot \frac{\varepsilon}{D} \cdot \mathrm{Re}\right)},\tag{5}$$

$$F''(\mathbf{x}\_i) = -\frac{172496738}{\ln(10) \cdot \left(9287 \cdot \mathbf{x}\_i + 1000 \cdot \frac{\varepsilon}{D} \cdot \mathbf{R}\varepsilon\right)^2},\tag{6}$$

For the first derivatives, [8,9,21–23] should be consulted.

For Example 1, in which *Re* = 2.3 <sup>×</sup> 105 and <sup>ε</sup> *<sup>D</sup>* = <sup>10</sup>−4, for initial starting point *x*<sup>0</sup> = √<sup>1</sup> *f*0 <sup>=</sup> <sup>ε</sup> *D* 3.71 + 1574003 ≈ 1574003, the final balanced solution will be achieved after 27 iterations; *x*<sup>27</sup> = √<sup>1</sup> *f*27 = 7.89313398624169 → *f*<sup>27</sup> = 0.016050961. However, for the initial starting point *x*<sup>0</sup> = √<sup>1</sup> *f*0 <sup>=</sup> <sup>ε</sup> *D* 3.71 − 2.4 ≈ −2.4, the procedure in MS Excel will report a division by zero error in the third step, so it will be a good opportunity for students to reexamine the importance of the appropriate chosen starting point. Thus, *x*<sup>0</sup> = √<sup>1</sup> *f*0 , which in this case can be selected from −1.172 to 1,583,778. Of course, negative starting points do not represent physical meaning. A similar exercise should be done for Example 2, in which *Re* = 4.6 <sup>×</sup> <sup>10</sup><sup>7</sup> and <sup>ε</sup> *<sup>D</sup>* = 0.037.

For this task, at least 90 min should be assigned, as multi-point iterative methods are more complex than two-point methods.

#### *2.2. Special Iterative Methods*

While simple operations, such as adding, subtracting, multiplication, and division, require minimal system resources of computers, more complex operations, such as logarithmic, non-integer power, and similar function calls, require multiple floating point operations to be executed in the Central Processing Unit (CPU), which represents a time- and energy-consuming task [15,25–28]. For example, students should understand that the execution of one non-integer power operation in a computer goes through execution of one logarithmic and one exponential function. Some numbers are too big or too small to be placed in registers of computers, which have 32 or 64 bits [37,38].

#### 2.2.1. Padé Approximant

As noted, one iterative cycle involves the evaluation of one logarithmic function when solving the Colebrook equation. The Colebrook's law is in its nature logarithmic, but here, students will learn how to avoid multiple evaluations of the logarithmic function by replacing it with a Padé approximant in all iterations, with the exception of the first one [30]. A Padé approximant is the "best" approximation of the given function by a rational function of given order. For example, log10(95) can be accurately approximated by its Padé approximant at the expansion point *x* = 1, keeping in mind that log10(100) <sup>=</sup> 2 and using the fact that log10(95) <sup>=</sup> log10- 100 100 95 <sup>=</sup> log10(100) <sup>−</sup> log10 100 95 where

100 <sup>95</sup> ≈ 1.0526 is close to 1. The Padé approximant will be very accurate in this case, because 1.0526 is very close to the expansion point *x* = 1.

Evaluation of the decadic logarithm can be expressed by the natural logarithm log10(*z*) <sup>=</sup> ln(*z*) ln(10), where *ln*(10) ≈ 2.302585093 is constant. The first task for students would be to make research related to Padé approximants. For students with interest in programming, the additional task would be to develop a code in appropriate software, such as Matlab. Padé approximants of different orders can be used for approximation of *ln*(*z*), for arguments *z* close to 1, i.e., *z* ≈ 1. As the expansion point *z*<sup>0</sup> = 1 is a root of *ln*(*z*), the accuracy of the Padé approximant decreases. Setting the "OrderMode" option in the Matlab Padé command to relative compensates for the decreased accuracy. Thus, here, the Pad<sup>é</sup> approximant of the (*m*,*n*) order uses the form *ln*(*z*) <sup>≈</sup> (*z*−*z*0)·(α0+α1·(*z*−*z*0)+...+α*m*(*z*−*z*0) *m*) 1+β1·(*z*−*z*0)+...+β*n*(*z*−*z*0) *<sup>n</sup>* , where α and β are coefficients (the coefficients of the polynomials need not be rational numbers). The Horner algorithm transforms polynomials into a computationally efficient form [30]. One possible procedure for solving the Colebrook equation using Padé approximants and the simple fixed-point iterative method is shown in Figure 2. Detailed step-by-step calculation using a table can be found in [30].

**Figure 2.** A possible procedure for solving the Colebrook equation using Padé approximants and the simple fixed-point iterative method.

Students interested in programming can use tic and toc functions ofMatlab and GNU/Octave [49,50] to measure the speed of calculations using the classical fixed-point method vs. the fixed-point method combined with Padé approximants (let us remind readers that GNU/Octave software is available for free [49]). Students can combine the Newton–Raphson method with Padé approximants and check the accuracy of Padé approximants of different orders [30].

For a general understanding of the concept of Padé approximants, 90 min is required (45 min for an introduction and a further 45 min for numerical tests). Then, that additional 90 min would be required for the application of Padé polynomials for solving the Colebrook equation (again, 45 min for an introduction and 45 min for numerical examples). Students with an interest in programming should have an additional 45-min exercise (a presentation of complexity vs. accuracy of Padé approximants).

#### 2.2.2. Approximations by Multi-Point Methods with Internal Cycles

Some iterative methods belong to the group of multi-point methods [8] that are very powerful and that can reach the accurate solution of the Colebrook equation in the first iteration. An example of an approximation of the Colebrook equation based on internal cycles is given by Serghides [51]. Serghides' methods belong to Steffensen types of methods, which do not require derivatives (other Steffensen types of methods can be used in addition [52,53]).

In order to develop an approximation with the internal cycle, students can choose one initial fixed point *<sup>f</sup>*0, for example *<sup>f</sup>*<sup>0</sup> <sup>=</sup> *fmin*<sup>+</sup> *fmax* <sup>2</sup> <sup>=</sup> 0.0064+0.077 <sup>2</sup> = 0.0417, and then calculate √1 *f*0 <sup>≈</sup> 4.9, which will give <sup>√</sup><sup>1</sup> *<sup>f</sup>* <sup>=</sup> <sup>−</sup>2· log10- 2.51·4.9 *Re* <sup>+</sup> <sup>ε</sup> *D* 3.71 and after one additional cycle √1 *<sup>f</sup>* ≈ −2· log10- ε *D* 3.71 <sup>−</sup> 5.02 *Re* · log10- 12.3 *Re* <sup>+</sup> <sup>ε</sup> *D* 3.71 . Using this exercise, students should see how important constant 12.3 is. Students should evaluate the relative error over the domain of applicability of the Colebrook's equation, then replace the constant 12.3 with a more appropriate value. Students can compare their chosen constant with values shown by Shacham [54] and by Zigrang and Sylvester [24]. For Example 1, in which *Re* = 2.3·105 and <sup>ε</sup> *<sup>D</sup>* = <sup>10</sup>−4, the approximate value is √1 *<sup>f</sup>* ≈ −2· log10 <sup>10</sup>−<sup>4</sup> 3.71 <sup>−</sup> 5.02 2.3·<sup>105</sup> · log10 12.3 2.3·10<sup>5</sup> <sup>+</sup> <sup>10</sup>−<sup>4</sup> 3.71 . Students should calculate this value and evaluate the introduced error compared with the accurate iterative solution [20]. An additional task could be to find additional approximations with internal cycles, such as Romeo et al. [55] or Offor and Alabi [56] and to discuss what additional strategies were used in these cases.

For this task, 45 min is required.

#### *2.3. Special Functions*

The Colebrook equation is implicit with respect to the unknown flow friction factor. However, the Colebrook equation can be rearranged in the explicit form using the Lambert W function, where further, this special function can be evaluated only approximately [5]. Some numerical constraints in using the Lambert W function, which can occur in our case, will be detected [37] and a way to mitigate this inconvenience will be shown [27–29].

For tasks related to special functions, students need 45 min for a general introduction and an additional 90 min to work on examples.

#### 2.3.1. Lambert W Function

Interested readers may find more information on the Lambert W function through [5,6,57,58]. Then, students should understand why some functions, such as logarithmic, square root, or exponential, are available in pocket calculators and spreadsheet solvers, while the Lambert W function is not [6]. In this task, students should understand how computers use a series expansion to interpret functions. Moreover, they need to understand how to use Visual Basic for Applications (VBA) modules in order to introduce new functions in MS Excel [59] (related to this task, students interested in programming should take at least 90 min).

The Lambert W function is used in mathematics to solve equations in which the unknown appears both outside and inside an exponential function or a logarithm, such as 3*x* + 2 = ex or *x* = ln(4·x). According to Fukushima [60], the Lambert W function is defined as the solution of a transcendental equation with <sup>κ</sup> = *W*(κ)·*eW*(κ) . Such equations cannot be, except in special cases, solved explicitly in terms of the finite number of arithmetic operations, nor in terms of exponentials or logarithms. In the case of the Colebrook equation, students can consult scientific literature, while here we will show formulation from [28,29], Equation (7), while a similar formulation can be found in [27]:

$$\begin{array}{ll} \frac{1}{\sqrt{f}} = \frac{2}{\ln(10)} \cdot \left( \ln \left( \frac{R\varepsilon}{2.51} \cdot \frac{\ln(10)}{2} \right) + \mathcal{W}(\varepsilon^x) - x \right) \\\ x = \ln \left( \frac{R\varepsilon}{2.51} \cdot \frac{\ln(10)}{2} \right) + \frac{R\varepsilon \cdot \frac{\varepsilon}{D}}{2.51 \cdot 3.71} \cdot \frac{\ln(10)}{2} \end{array} \tag{7}$$

Our carefully chosen examples, *Re* = 2.3 <sup>×</sup> 105 and <sup>ε</sup> *<sup>D</sup>* <sup>=</sup> <sup>10</sup>−4, (2) *Re* <sup>=</sup> 4.6 <sup>×</sup> 107 and <sup>ε</sup> *<sup>D</sup>* = 0.037, respectively, are chosen to show that using floating-point arithmetic, Equation (7), gives the accurate solution for Example 1, while for Example 2, an arithmetic overflow error will occur. It is because *W*(*ex*) becomes too large for ordinary computer registers.

#### 2.3.2. Wright ω Function and Related Approximations

As shown in Biberg [27] and Brki´c and Praks [28,29], the problem with too large numbers in computer registers can be solved by introducing a cognate of the Lambert W function—the Wright ω function. The Wright ω function *y* = ω(*x*) = *W*(*ex*) is a solution to the equation *y* + *ln*(*y*) = *x*. Here, students should reformulate Equation (7) through the Wright ω function [50].

For Example 2, which uses *Re* = 4.6 <sup>×</sup> 107 and <sup>ε</sup> *<sup>D</sup>* = 0.037, *x* is approximately 210441.7. In Matlab, the command *y* = *wrightOmega*(210441.7) returns the value *y* = 210429.44.

In Matlab, but also in the GNU/Octave software [49], students can use the free WrightOmegaq library [50]. In this case, the command is *y* = *wrightOmegaq*(210441.7).

Students can easily verify that the solution *y* satisfies the relation *y* + *ln*(*y*) = *x*. Let us reiterate that for Example 2, the value of *x* is approximately 210441.7.

In the following task, students should explore how to estimate the Wright ω function without Matlab or GNU/Octave. The simplification *W*(*ex*) − *x* ≈ *ln*(*x*)· 1 *<sup>x</sup>* − 1 gives an accurate approximation suitable for engineering practice—Equation (8) [28,29]:

$$\begin{array}{c} \frac{1}{\sqrt{f}} \approx 0.8686 \cdot \left[ B + \ln(B+A) \cdot \left( \frac{1}{B+A} - 1 \right) \right] \\\ A \approx \frac{Rc \cdot \frac{1}{D}}{8.0078} \\\ B \approx \ln \left( \frac{Rc}{2.18} \right) \approx \ln(Re) - 0.7794 \end{array} \tag{8}$$

where <sup>2</sup>·2.51 *ln*(10) ≈ 2.18 and 2.18·3.71 ≈ 8.0878.

Students should evaluate the relative error of the approximation, Equation (8), and develop a more accurate approximation using the enhanced model *W*(*ex*) − *x* ≈ *c* + *ln*(*x*)· 1 *<sup>x</sup>* − 1 , where *c* is a constant. Of course, if *c* = 0, the enhanced model is equivalent to the model (8). For the general case, the constant *c* can be found by the MS Excel Solver, in order to minimize the relative error of the model with the iterative solution of the Colebrook equation provided by MS Excel. Students should do this optimization task with randomly selected input values of the Colebrook model *Re*, <sup>ε</sup> *D* , which are sampled from the domain of applicability of the Colebrook equation. Finally, students can compare the relative error of this enhanced model with the relative error of the original model.

#### 2.3.3. Tania Function and Related Approximations

Students should examine the Tania function, which is defined as *Tania*(*x*) = ω(*x* + 1), *Tania*(*x*) = *W*(*ex*), all with reference to [61] (the Euler T function should be an appropriate task to be further examined [39]).

For this task, 45–90 min should be assigned, during which students need to try to develop ability accepting new functions with new concepts [62].

#### **3. Conclusions**

Here we present how the relatively simple Colebrook equation, widely known by all students of hydraulics, can be used for introducing in curricula various iterative methods, special functions such as Lambert W, polynomials such as Padé approximants, optimization methods, etc.

The proposed exercises are suitable for students of hydraulics (civil, mechanical, chemical, or petroleum engineering curriculum, including the exploitation of oil and gas), but can be used as practical examples in teaching of numerical methods for students of mathematics. The proposed exercises can be extended for the teaching of pipe network analysis [31–36]. These methods can be combined with recent experimental works [63]. Furthermore, instead of pipes, a modified version of the Colebrook equation can be used for fuel cells [64,65].

**Author Contributions:** The paper is a product of the joint efforts of the authors, who worked together on models of natural gas distribution networks. P.P. has a scientific background in applied mathematics and programming, while D.B. has a background in control and applied computing in mechanical and petroleum engineering. They used the experience gained from their previous papers related to the Colebrook equation to develop proposed teaching exercises. D.B. wrote a draft of this article and prepared figures.

**Funding:** This work has been partially funded by the Technology Agency of the Czech Republic, partially by the National Centre for Energy TN01000007 and partially by the project "Energy System for Grids" TK02030039. Work of D.B. has been also supported by the Ministry of Education, Youth and Sports from the National Programme of Sustainability (NPS II) project "IT4Innovations excellence in science - LQ1602" and by the Ministry of Education, Science, and Technological Development of the Republic of Serbia through the project "Development of new information and communication technologies, based on advanced mathematical methods, with applications in medicine, telecommunications, power systems, protection of national heritage and education" iii44006.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Teach Second Law of Thermodynamics via Analysis of Flow through Packed Beds and Consolidated Porous Media**

#### **Rajinder Pal**

Department of Chemical Engineering, University of Waterloo, Waterloo, ON N2L 3G1, Canada; rpal@uwaterloo.ca; Tel.: +1-519-888-4567 (ext. 32985)

Received: 12 June 2019; Accepted: 24 June 2019; Published: 27 June 2019

**Abstract:** The second law of thermodynamics is indispensable in engineering applications. It allows us to determine if a given process is feasible or not, and if the given process is feasible, how efficient or inefficient is the process. Thus, the second law plays a key role in the design and operation of engineering processes, such as steam power plants and refrigeration processes. Nevertheless students often find the second law and its applications most difficult to comprehend. The second law revolves around the concepts of entropy and entropy generation. The feasibility of a process and its efficiency are directly related to entropy generation in the process. As entropy generation occurs in all flow processes due to friction in fluids, fluid mechanics can be used as a tool to teach the second law of thermodynamics and related concepts to students. In this article, flow through packed beds and consolidated porous media is analyzed in terms of entropy generation. The link between entropy generation and mechanical energy dissipation is established in such flows in terms of the directly measurable quantities such as pressure drop. Equations are developed to predict the entropy generation rates in terms of superficial fluid velocity, porous medium characteristics, and fluid properties. The predictions of the proposed equations are presented and discussed. Factors affecting the rate of entropy generation in flow through packed beds and consolidated porous media are identified and explained.

**Keywords:** undergraduate education; applications of fluids; fluid mechanics; packed bed; porous media; non-equilibrium thermodynamics; entropy generation; pressure loss; Ergun equation; Forchheimer equation

#### **1. Introduction**

Thermodynamics is a difficult subject to learn and teach. It is considered to be one of the most abstract disciplines of the physical sciences [1]. Students all over the world face difficulties in learning thermodynamics. More specifically, it is the second law of thermodynamics dealing with entropy and entropy production that is difficult for students to fully comprehend. The second law of thermodynamics simply states that all real (irreversible) processes are accompanied by the production of entropy in the universe. However, the students of thermodynamics often find it difficult to visualize entropy generation in real processes at a mechanistic level. The quantification of entropy generation in real processes is an equally problematic issue for students. What causes entropy generation and how do we quantify entropy generation? These are some of the fundamental questions faced by students. Moreover, there are no instruments which can be used to directly measure entropy and entropy generation in real processes.

Entropy generation occurs in flow of all real (viscous) fluids [2]. When a viscous fluid is forced to flow through any geometry, such as a pipe, viscous stresses and velocity gradients are established. For fluid to remain in motion, work has to be done against the viscous stresses which oppose its motion.

Consequently, part of the mechanical energy of fluid is dissipated into frictional heat (internal energy) during motion. The dissipation of highly ordered mechanical energy into disorderly internal energy is reflected in entropy generation. Thus, the analysis of fluid mechanics problems and measurement of the appropriate flow variables such as pressure loss could be used as a tool to demonstrate and quantify entropy generation in real processes.

The main objectives of this article are: (a) to analyze the flow of viscous fluids through packed beds of discrete particles and through consolidated porous media in terms of entropy generation; (b) to quantify entropy generation in flow through packed beds and consolidated porous media in terms of directly measureable quantities such as pressure loss as a function of flow rate of fluid; and (c) to discuss various factors which affect the rate of entropy generation in flow through packed beds and consolidated porous media.

#### **2. Brief Review of the Second Law of Thermodynamics**

The second law of thermodynamics could be stated in several different but equivalent ways. The classical statements of the second law are the Kelvin–Planck statement and the Clausius statement.

The Kelvin–Planck statement says that "It is impossible for a system operating in a cycle and connected to a single heat reservoir to produce a positive amount of work in the surroundings [3]". According to the Kelvin–Planck statement, it is impossible to build a heat engine shown schematically in Figure 1 where heat absorbed from a heat reservoir is completely converted into work without altering the properties of the system. Mathematically, the Kelvin–Planck statement can be expressed as [4]:

$$\mathcal{W}\_{cycle} \le 0 \text{ (single heat reservoir)}\tag{1}$$

where the system communicates thermally only with a single heat reservoir. The sign convention for work (*W*) used in this article is *W* = positive, if it is produced by the system (flows out of the system to surroundings) and *W* = negative, if it is absorbed by the system (flows into the system from surroundings). Thus, no cyclic process is possible where *Wcycle* > 0 using a single heat reservoir. The equality in Equation (1) is valid for a reversible process and the inequality is valid for an irreversible process. For a reversible process, no work is produced or destroyed, that is, *Wcycle* = 0 whereas work is destroyed when the process is irreversible, that is, *Wcycle* < 0 (assuming that the system communicates thermally only with a single heat reservoir).

**Figure 1.** Impossible heat engine.

The Clausius statement of the second law [3] says that "It is impossible for a system operating in a cycle to have its *sole e*ff*ect* the transfer of heat from a low temperature heat reservoir to a high temperature heat reservoir". According to the Clausius statement, it is impossible to construct a device based on the scheme shown schematically in Figure 2 where the sole result of the process is the transfer of heat from a cooler body at low temperature (*TL*) to a hotter body at high temperature (*TH*).

**Figure 2.** Impossible scheme.

Another powerful statement of the second law of thermodynamics is that all irreversible (real) processes are accompanied by entropy generation in the universe [5], that is:

$$S\_{G,univ} = \Delta S\_{sys} + \Delta S\_{surr} \ge 0 \tag{2}$$

where *SG*,*univ* is the total amount of entropy generated in the universe (system + surroundings), Δ*Ssys* is the entropy change of the system and Δ*Ssurr* is the entropy change of the surroundings. The equality in Equation (2) is valid for a reversible process and the inequality is valid for an irreversible process. According to this statement of the second law, no process is possible for which *SG*,*univ* < 0.

Entropy is a measure of disorderliness of a system. According to the Boltzmann entropy equation, the entropy of a system can be expressed as [3]:

$$S = k\_{\rm B} \ln \Omega \tag{3}$$

where Ω is the number of possible configurations of the system and *kB* is the Boltzmann constant. The larger the number of possible configurations of the system, greater is the disorderliness of the system and higher is the entropy. Therefore, we can interpret the second law of thermodynamics in

yet another way, that is, "Only those processes are possible processes which lead to an increase in the disorderliness of the universe".

The scheme shown in Figure 1 is impossible as it decreases the disorder of the universe, that is, *SG*,*univ* < 0. Here Δ*Ssys* = 0 but Δ*Ssurr* < 0 as the surroundings heat reservoir loses heat. For this scheme, Equation (2) gives:

$$S\_{G, \text{univ}} = -\frac{Q\_{\text{cyclr}}}{T} \ge 0 \tag{4}$$

where *T* is the absolute temperature of the heat reservoir and *Q* is the heat transferred. As *Qcycle* = *Wcycle* from the first law of thermodynamics, Equation (4) reduces to Equation (1), that is, *Wcycle* ≤ 0 when there is only one heat reservoir involved. Similarly, the scheme shown in Figure 2 is impossible as it decreases the disorder of the universe:

$$
\Delta S\_{\text{G,univ}} = \Delta S\_{\text{cold}-body} + \Delta S\_{\text{hot}-body} < 0 \tag{5}
$$

Note that Δ*Scold*−*body* < 0 as it loses heat whereas Δ*Shot*−*body* > 0 as it gains heat. However, due to different temperatures of the cold and hot bodies, Δ*Scold*−*body* > Δ*Shot*−*body* .

Thus, the Kelvin–Planck and Clausius statements of the second law of thermodynamics are special cases of the statement of the second law expressed in the form of Equation (2).

For a flow process (see Figure 3), the second law of thermodynamics can be written as [2]:

$$\dot{S}\_{\text{G,min}} = \dot{S}\_{\text{G,C}} + \dot{S}\_{\text{G,Surr}} = \iiint \rho s(\hat{n} \cdot \vec{V}) dA + \iiint \frac{\partial(\rho s)}{\partial t} d\theta - \sum \frac{\dot{Q}\_i}{T\_i} \ge 0 \tag{6}$$

where . *SG* is the rate of entropy generation, *s* is the entropy per unit mass of fluid, ρ is the fluid density, *<sup>n</sup>*<sup>ˆ</sup> is unit outward normal to the control surface, <sup>→</sup> *V* is fluid velocity vector, *A* is the control surface area, *<sup>t</sup>* is time, <sup>ϑ</sup> is the volume of the control volume, . *Qi* is the rate of heat transfer to control volume from *i*th heat reservoir at an absolute temperature of *Ti*, the subscripts *CV* and *Surr* refer to control volume and surroundings, respectively. As noted earlier, the equality in Equation (6) is valid for a reversible (frictionless) process and the inequality is valid for an irreversible process. The surface integral <sup>ρ</sup>*s*(*n*ˆ· → *V*)*dA* is the net outward flow of entropy (associated with mass) across the entire control surface. The volume integral <sup>∂</sup>(ρ*s*) <sup>∂</sup>*<sup>t</sup> d*ϑ is the rate of accumulation of entropy within the entire control volume (assumed to be fixed and non-deforming). For a control volume with one inlet and one outlet, Equation (6), under steady state condition, reduces to:

$$\dot{S}\_{\text{G,universal}} = \dot{m}(\Delta s) - \sum \frac{\dot{Q}\_i}{T\_i} \ge 0 \tag{7}$$

where . *m* is the mass flow rate.

The quantification of entropy generation in real processes is important from a practical point of view as entropy generation is directly related to the efficiency of the process. Higher the rate of entropy generation in a process, lower is the thermodynamic efficiency of the process. According to the Gouy–Stodola theorem [6] of thermodynamics, the loss of power or work potential in a real process, due to irreversibilities in the process, is directly proportional to the total rate of entropy generation. Thus: .

$$
\dot{\mathcal{W}}\_{\text{lost}} \propto \dot{\mathcal{S}}\_{\text{G},\text{univ}} \tag{8}
$$

where . *Wlost* is the rate of work lost (wasted) as a result of irreversibilities in the process.

As an example of a flow process, consider flow through a control volume shown in Figure 3. The first law of thermodynamics for open systems under steady state condition gives:

$$
\dot{m} \left[ \Delta \text{h} + \Delta \left( V^2 / 2 \right) + \text{g} \, \Delta z \right] = \dot{Q} - \dot{W}\_{\text{sh}} \tag{9}
$$

where *h* is the specific enthalpy of fluid, *V* is the fluid velocity, *g* is the acceleration due to gravity, *z* is the elevation, . *<sup>Q</sup>* is the rate of heat transfer, and . *Wsh* is the rate of shaft work. Neglecting kinetic and potential energy changes, Equation (9) simplifies to:

$$
\dot{\mathcal{W}}\_{sh} = \dot{\mathcal{Q}} - \dot{m}\Delta h = \dot{m}(\mathcal{Q} - \Delta h) \tag{10}
$$

where *Q* is the heat transfer per unit mass of fluid. The second law, Equation (7), can be written as:

$$
\dot{S}\_{\text{G,univ}} = \dot{m}(\Delta s) - \frac{\dot{Q}}{T\_o} \ge 0 \tag{11}
$$

where *To* is the absolute temperature of the heat reservoir (surroundings). For the flow process to be reversible: .

$$\dot{S}\_{\text{G,univ}} = \dot{m} \left[ (\Delta s) - \frac{Q\_{\text{rev}}}{T\_o} \right] = 0 \tag{12}$$

$$Q\_{rev} = T\_o(\Delta \text{s}) \tag{13}$$

where *Qrev* is the heat transfer per unit mass of fluid for a reversible process. From Equations (10) and (13): .

$$
\dot{\mathcal{W}}\_{sl, rev} = \dot{m} \left[ T\_o(\Delta \mathbf{s}) - \Delta \mathbf{h} \right] \tag{14}
$$

The lost work, . *Wlost*, is defined as:

$$
\dot{\mathcal{W}}\_{\rm lost} = \dot{\mathcal{W}}\_{\rm sh, rev} - \dot{\mathcal{W}}\_{\rm sh} \tag{15}
$$

From Equations (10), (14), and (15), we get:

$$
\dot{W}\_{\text{lost}} = \dot{m} [T\_o(\Delta \text{s}) - Q] = \dot{m} T\_o \left[ (\Delta \text{s}) - \frac{Q}{T\_o} \right] \tag{16}
$$

Using Equation (11), Equation (16) can be re-written as:

$$
\dot{\mathcal{W}}\_{\text{lost}} = T\_o \dot{\mathcal{S}}\_{\text{G,\omega miv}} \ge 0 \tag{17}
$$

Equation (17) is the Gouy–Stodola theorem [6] of thermodynamics. Thus, the rate of work lost due to irreversibilities is directly proportional to the rate of entropy generation. The thermodynamic efficiency of a flow process can be defined as:

$$\varepsilon = \frac{\dot{W}\_{sh}}{\dot{W}\_{sh, rev}} = 1 - \frac{\dot{W}\_{lost}}{\dot{W}\_{sh, rev}} \le 1.0 \tag{18}$$

In a reversible process, . *Wlost* <sup>=</sup> 0 and <sup>=</sup> 1. In any irreversible process, . *Wlost* > 0 and < 1. If no work is produced in the process, that is, actual . *Wsh* is zero, then all the work potential is lost due to irreversibilities in the process and consequently, = 0.

**Figure 3.** Typical flow process.

#### **3. Flow through Unconsolidated and Consolidated Porous Media**

Porous medium is a composite material in that it consists of two phases, namely pores (voids, free space pockets) and solid-phase. The pores may be occupied by a fluid (gas, oil, water, etc.). A large variety of natural and synthetic materials are porous in nature. Examples include: underground oil reservoirs, ceramics, solid foams, sand filters, wood, and packed beds of particles used widely in chemical engineering applications. The pores of a porous medium usually form a three-dimensional inter-connected network and, therefore, fluids can flow through the porous medium. If all pores of a porous medium are inter-connected, then the porosity (ε) of the porous medium is simply the fraction of the total volume of the medium that is occupied by the pores. Thus, the fraction of the total volume that is occupied by the solid phase is (1 − ε). When some pores are isolated or disconnected or have dead ends, then the effective porosity, defined as the ratio of connected void volume to total volume of the medium, is lower than the total porosity.

Porous media could be classified as consolidated or unconsolidated. In consolidated porous medium, the solid phase is basically a single piece of material or the grains of the solid phase are cemented or fused together to form a single piece of solid phase. In unconsolidated porous media, on the other hand, the grains or the particles of the solid phase are not cemented together and, therefore, the porous medium is a multi-particle system like a packed bed of individual (un-cemented, un-glued) particles. Flow of single-phase Newtonian fluids (gas, water, oil) through packed beds and consolidated porous media is important from a practical point of view.

#### *3.1. Analysis of Flow through Packed Beds (Unconsolidated Porous Media)*

#### 3.1.1. Pressure-Loss in Flow through Packed Beds

Consider flow of a Newtonian fluid through a horizontal cylindrical packed bed of particles with bed diameter *Db* and bed length *Lb*. As flow through a packed bed is quite complex, a rigorous theoretical derivation of pressure drop-flow rate relationship is not possible. Only approximate models have been developed. In one approach, used widely to model flow through packed bed of particles, the packed bed is visualized as a bundle of identical capillary tubes [7–9]. In its simplest form, the capillary tube bundle model assumes that the capillary tubes are straight, cylindrical of constant cross-section (uniform radius), and parallel (all oriented in the same direction), as shown in Figure 4. The average velocity in any capillary tube is the same as that in the packed bed. The capillary diameter is equal to the hydraulic diameter of the bed, *dH*, defined as:

$$d\_{H} = 4 \left( \frac{volume\ of\ voids}{total\ wetted\ surface\ area} \right) = \left(\frac{4\varepsilon}{1-\varepsilon}\right) \left(\frac{V\_{p}}{S\_{p}}\right) \tag{19}$$

where ε is the void fraction or porosity, *Vp* is the volume of a single particle, and *Sp* is the surface area of a single particle. In writing the above expression for *dH*, it is assumed that all the particles are identical and that the wetted surface of the cylindrical container walls of the bed is negligible as compared with the total wetted surface area of the particles. For a bed of identical spherical particles, *Vp*/*Sp* = *dp*/6, where *dp* is the particle diameter. Thus, the hydraulic diameter becomes:

$$d\_{H} = \left(\frac{2}{3}\right) \left(\frac{\varepsilon}{1-\varepsilon}\right) d\_{p} \tag{20}$$

**Figure 4.** Path taken by fluid in a bed of particles and tube bundle model of packed bed.

The superficial velocity of fluid (*Vs*) is defined is:

$$V\_s = \frac{\dot{\mathfrak{g}}}{A\_b} \tag{21}$$

where . ϑ is the volumetric flow rate of fluid and *Ab* is the total cross-sectional area of the bed. Thus, *Vs* is the velocity of fluid in the bed if no particles were present in the bed. The average velocity *V* of fluid in the bed, also called interstitial velocity, is defined as:

$$
\overline{V} = \frac{\dot{\mathfrak{g}}}{A\_{flow}} \tag{22}
$$

where *Aflow* is the cross-sectional area of bed through which the fluid flows. The porosity of the bed ε is given as:

$$
\kappa = \frac{A\_{flow}}{A\_b} \tag{23}
$$

Therefore, the average or interstitial velocity of fluid in the bed is:

$$
\overline{V} = \frac{V\_s}{\varepsilon} \tag{24}
$$

This expression of average velocity does not consider the tortuous path taken by fluid in the bed. Due to tortuosity of the bed, the actual average velocity of fluid in the bed is larger than that given by Equation (24). The tortuosity τ of a bed is defined as:

$$
\pi = \frac{L\_{\prec}}{L\_{b}} \tag{25}
$$

where *Le* is average length of the tortuous path taken by the fluid and *Lb* is the straight length of the bed. To account for the tortuosity, the actual average velocity of the fluid in the bed can be expressed as follows:

$$
\overline{V} = \pi \frac{V\_s}{\varepsilon} \tag{26}
$$

In laminar flow of a Newtonian fluid through a cylindrical tube, the pressure gradient in the direction of flow is given as:

$$-\left(\frac{dP}{d\mathbf{x}}\right) = 32\left(\frac{\mu\overline{V}}{D^2}\right) \tag{27}$$

where μ is the fluid viscosity and *D* is the tube diameter. Replacing *D* with hydraulic diameter *dH* and the constant factor of 32 with *C*, the pressure-gradient in a capillary tube model of the packed-bed model can be expressed as:

$$-\left(\frac{dP}{d\mathbf{x}}\right) = \mathbf{C}\left(\frac{\mu\overline{V}}{d\_H^2}\right) \tag{28}$$

Note that flow in the bed is assumed to be laminar here. The constant *C* is expected to be larger than 32, as the path of fluid in the bed is not straight. The fluid follows a tortuous path in the bed and consequently, the pressure drop over a certain straight length is expected to be more than that observed over the same length if the fluid path was non-tortuous. Therefore, the constant *C* is taken as 32 τ where τ is the tortuosity. Upon substitution of the expressions for *dH* and *V*, and taking C = 32 τ, Equation (28) gives:

$$-\left(\frac{dP}{d\mathbf{x}}\right) = 72\,\,\tau^2 \left(\frac{\mu V\_s}{d\_p^2}\right) \frac{\left(1-\varepsilon\right)^2}{\varepsilon^3} \tag{29}$$

Equation (29) assumes that the cross-section of the representative flow passage (capillary tube) in the porous medium is circular and constant. In reality the fluid moves through converging-diverging flow passages of non-uniform cross-sections. In converging-diverging flows, fluid experiences stretching or extensional deformation. The elongation or stretching of fluid elements results in additional dissipation of energy and pressure drop. Thus, this equation needs to be modified further as follows:

$$-\left(\frac{dP}{d\mathbf{x}}\right) = 72\,\mathrm{\tau}^2\mathrm{K}\left(\frac{\mu V\_s}{d\_p^2}\right)\frac{\left(1-\varepsilon\right)^2}{\varepsilon^3} \tag{30}$$

where *K* is an empirical factor that takes into account the influence of non-constant cross-section of flow passage on pressure-gradient. The tortuosity <sup>τ</sup> is often taken to be <sup>√</sup> 2 for random packing of uniform spheres [9]. It is generally a function of the porosity of the porous medium [10]. For example, the following model, based on the Maxwell equation for electrical conductivity of composite, is often used to describe the relationship between tortuosity and porosity [10]:

$$
\pi = 1.5 - 0.5\varepsilon \tag{31}
$$

When <sup>τ</sup> <sup>=</sup> <sup>√</sup> 2 and *K* = 25/24, the following well-known Blake–Kozeny equation for laminar flow through randomly packed bed of uniform spheres is obtained [7]:

$$-\left(\frac{dP}{d\mathbf{x}}\right) = 150 \left(\frac{\mu V\_s}{d\_p^2}\right) \frac{\left(1 - \varepsilon\right)^2}{\varepsilon^3} \tag{32}$$

When <sup>τ</sup> <sup>=</sup> <sup>√</sup> 2 and *K* = 5/4, the following Carman–Kozeny equation, another well-known equation for laminar flow through randomly packed bed of uniform spheres, is obtained [11]:

$$-\left(\frac{dP}{d\mathbf{x}}\right) = 180 \left(\frac{\mu V\_s}{d\_p^2}\right) \frac{\left(1 - \varepsilon\right)^2}{\varepsilon^3} \tag{33}$$

Both Blake—Kozeny and Carman–Kozeny equations are popular in the literature although they use different values of the factor *K*. The observed difference in *K* values is probably related to differences in particle-shape, surface roughness of particles, and porosity of bed. Note that these equations are restricted to only laminar flow through bed of nearly uniform spheres. The packed-bed Reynolds number *Reb*, defined below, should be less than 10 for flow to be laminar in the bed.

$$R\varepsilon\_{b} = \left(\frac{1}{1-\varepsilon}\right) \frac{\rho d\_{p} V\_{s}}{\mu} \tag{34}$$

For turbulent flow through packed beds (*Reb* > 1000), the following Burke–Plummer equation is often used [12]:

$$1 - \left(\frac{dP}{dx}\right) = 1.75 \left(\frac{\rho V\_s^2}{d\_p}\right) \frac{(1 - \varepsilon)}{\varepsilon^3} \tag{35}$$

In the transition region, the following equation, obtained by superposition of Blake–Kozeny and Burke–Plummer equations, proposed by Ergun [13] is widely used:

$$-\left(\frac{dP}{d\mathbf{x}}\right) = 150 \left(\frac{\mu V\_s}{d\_p^2}\right) \frac{\left(1 - \varepsilon\right)^2}{\varepsilon^3} + 1.75 \left(\frac{\rho V\_s^2}{d\_p}\right) \frac{\left(1 - \varepsilon\right)}{\varepsilon^3} \tag{36}$$

As the Ergun equation, Equation (36), is obtained by the superposition of laminar and turbulent expressions, it is valid over the full range of the packed-bed Reynolds number *Reb*. The Ergun equation could also be re-cast in the following form [14]:

$$f\_b = \frac{150}{\text{Re}\_b} + 1.75\tag{37}$$

where *fb* is the packed-bed friction factor defined as:

$$f\_b = \frac{-\left(\frac{dP}{dx}\right)}{\left(\frac{\rho V\_s^2}{d\_p}\right)\frac{(1-\varepsilon)}{\varepsilon^3}}\tag{38}$$

The Ergun equation, Equation (36) or (37), is used extensively in the literature to describe pressure loss in packed beds. However, the following points should be kept in mind when using the Ergun equation: (a) it is applicable to unconsolidated beds of nearly uniform-size spherical particles with appreciable roughness. For smooth spheres, it tends to over-predict *fb* in the high *Reb* range (*Reb* > 700) [14–16]; (b) if the bed particles are non-uniform in sizes, the Sauter mean diameter of particles should be used in the application of the Ergun equation; (c) when particle shape deviates significantly from a sphere, the Ergun equation tends to under-predict *fb* [17]; and (d) wall effects can be important when *Db*/*dp* < 10. When wall effects are important, the experimental data show deviation from the predictions of the Ergun equation [18].

Figure 5 shows the prediction of the packed bed friction factor as a function of packed bed Reynolds number using the Ergun equation, Equation (37). In the limits of low *Reb* and high *Reb*, the Ergun equation predictions overlap with the predictions of the Blake–Kozeny equation (low *Reb*) and the Burke–Plummer equation (high *Reb*), as expected.

**Figure 5.** Prediction of packed bed friction factor as a function of packed bed Reynolds number. The data points are generated using the Ergun equation (Equation (37)).

#### 3.1.2. Entropy Generation in Flow through Packed Beds

Consider steady flow of an incompressible fluid through a bed of particles (see Figure 6). We can apply the following macroscopic mechanical energy balance between the inlet and outlet of the bed [2]:

$$
\dot{\mathcal{W}}\_{\rm sh} + \dot{F}\_l + \dot{m} \left[ \Delta \varrho + \Delta (\rm KE) + \frac{1}{\rho} (\Delta P) \right] = 0 \tag{39}
$$

where . *Fl* is the rate of mechanical energy dissipation due to friction in fluid, Δϕ is the potential energy change per unit mass of fluid, Δ(*KE*) is the kinetic energy change per unit mass of fluid, and Δ*P* is the pressure change of the fluid. Neglecting potential and kinetic energy changes and taking . *Wsh* = 0, the mechanical energy balance gives:

$$
\dot{F}\_I = -\frac{\dot{m}}{\rho}(\Delta P) \tag{40}
$$

The rate of mechanical energy dissipation per unit length of packed bed ( . *F <sup>l</sup>*) can be expressed as:

$$
\dot{F}\_I' = \frac{\dot{m}}{\rho} \left( -\frac{dP}{dx} \right) \tag{41}
$$

From the second law of thermodynamics, Equation (11), we get:

$$
\dot{S}\_{G,\mu\dot{m}\dot{v}} = \dot{m}(\Delta s) > 0\tag{42}
$$

Note that we are assuming flow to be adiabatic with negligible heat transfer. There is no entropy generation in the surroundings. All the entropy is generated within the fluid inside the packed bed and the rate of entropy generation is the net rate of increase in entropy of the flowing stream.

**Figure 6.** Flow through a packed bed

We can now relate entropy change of the fluid stream to pressure change. For pure substances, the relationship between entropy and other state variables is given as [5]:

$$Tds = dh - \left(dP/\rho\right)\tag{43}$$

where *T* is the absolute temperature. From the first law of thermodynamics, Equation (9), the enthalpy change is zero in the absence of heat transfer and shaft work for steady flow in a horizontal bed. Consequently, Equation (43) reduces to:

$$Tds = -\left(dP/\rho\right)\tag{44}$$

Assuming incompressible flow and constant temperature, Equation (44) upon integration gives:

$$
\Delta s = -\frac{\Delta P}{\rho T} \tag{45}
$$

From Equations (42) and (45), it follows that:

$$
\dot{S}\_G = \frac{\dot{m}}{\rho} \left( -\frac{\Delta P}{T} \right) > 0 \tag{46}
$$

The subscript "*univ*" has been removed from . *SG*,*univ* as . *SG*,*univ* is simply the rate of entropy generation within the fluid inside a packed bed. From Equations (41) and (46), we can also express the rate of entropy generation in a packed bed on a unit volume basis as:

$$
\dot{S}\_{G}^{\prime\prime\prime} = \frac{\dot{F}\_{l}^{\prime}}{TA\_{b}} = \frac{V\_{s}}{T} \left(-\frac{dP}{dx}\right) \tag{47}
$$

where *Ab* is the total cross-sectional area of the bed and . *S <sup>G</sup>* is the rate of entropy generation per unit volume of the bed. From Equations (38) and (47), we get:

$$\dot{S}\_{G}^{\prime\prime\prime} = \left(\frac{1}{T}\right) \left(\frac{1-\varepsilon}{\varepsilon^3}\right) \left(\frac{\rho V\_s^3}{d\_p}\right) f\_b \tag{48}$$

The packed bed friction factor *fb* for laminar flow (*Reb* < 10) is given as:

$$f\_b = \frac{150}{Rc\_b} \tag{49}$$

Consequently, Equation (48) reduces to:

$$\dot{S}\_{G}^{\prime\prime} = \left(\frac{1}{T}\right) \left(\frac{150\left(1-\varepsilon\right)^{2}}{d\_{p}^{2}\varepsilon^{3}}\right) \mu V\_{s}^{2} \tag{50}$$

Thus, entropy generation rate per unit volume of bed in steady laminar flow of a Newtonian fluid is directly proportional to the fluid viscosity and to the square of the superficial fluid velocity in the bed. The entropy generation rate also depends on the particle diameter, the bed porosity ε, and the temperature.

The packed bed friction factor in turbulent flow (*Reb* > 1000) is given as [7–9]:

$$f\_b = 1.75\tag{51}$$

Consequently, Equation (48) reduces to:

$$\dot{S}\_{G}^{\prime\prime} = \left(\frac{1}{T}\right) \left(\frac{1.75(1-\varepsilon)}{d\_P \varepsilon^3}\right) \left(\rho V\_s^3\right) \tag{52}$$

Thus, entropy generation rate per unit volume of bed in steady turbulent flow of a Newtonian fluid is independent of the fluid viscosity and is directly proportional to the cube of the superficial fluid velocity in the bed. The entropy generation rate also depends on the fluid density, the particle diameter, the bed porosity ε, and the temperature.

To cover the full range of packed bed Reynolds number *Reb*, we substitute the Ergun equation, Equation (37), into Equation (48) to obtain the following equation valid over the full range of packed bed Reynolds number *Reb*:

$$\dot{S}\_{G}^{\prime\prime} = \frac{1}{T} \Biggl[ \left( \frac{150(1-\varepsilon)^2}{d\_p^2 \varepsilon^3} \right) \mu V\_s^2 + \left( \frac{1.75(1-\varepsilon)}{d\_p \varepsilon^3} \right) \left( \rho V\_s^3 \right) \right] \tag{53}$$

Note that Equation (53) is simply the superposition of Equations (50) and (52).

#### *3.2. Analysis of Flow through Consolidated Porous Media*

#### 3.2.1. Pressure-Loss in Flow through Consolidated Porous Media

Laminar flow in consolidated porous media (see Figure 7) is usually described by Darcy's law given as follows [19]:

$$-\frac{dP}{dx} = \frac{\mu}{k}V\_s\tag{54}$$

where *k* is the permeability of the porous medium. Permeability is a measure of the ability of porous medium to conduct fluid. Higher *k* means lower resistance to flow and consequently, larger flow rate for the same pressure gradient. *k* has the dimensions of (length)2. It is usually expressed in darcies. For example, 1 darcy <sup>=</sup> 1 (cm/s). centipoise/(atm/cm) <sup>=</sup> 9.87 <sup>×</sup> <sup>10</sup>−<sup>13</sup> <sup>m</sup>2.

**Figure 7.** Consolidated porous medium (adapted from www.fotosearch.com.au).

Darcy's law is equally applicable to unconsolidated porous medium such as packed bed of free (un-cemented/un-fused) particles. Thus, the Blake–Kozeny equation or the Carman–Kozeny equation are special cases of the Darcy law. Upon comparison of the Darcy law with Blake–Kozeny/ Carman–Kozeny equations, the permeability of a packed bed of uniform size spherical particles can be expressed as follows:

$$k = \frac{d\_p^2 \varepsilon^3}{150(1 - \varepsilon)^2} \tag{55}$$

$$k = \frac{d\_p^2 \varepsilon^3}{180(1 - \varepsilon)^2} \tag{56}$$

Equation (55) gives permeability based on the Blake–Kozeny equation whereas Equation (56) gives permeability based on the Carman–Kozeny equation.

Darcy's law is restricted to flows where viscous forces dominate over the inertial forces. At high flow rates, the flow in the porous medium becomes turbulent and Darcy's law is no longer valid. The pressure-gradient in the turbulent regime is much higher compared with that in the laminar regime at the same superficial velocity *Vs*. In order to extend the Darcy law to turbulent regime, Forchheimer [19] modified the Darcy law by adding a quadratic term as follows:

$$-\frac{dP}{d\mathbf{x}} = \frac{\mu}{k}V\_s + \beta\rho V\_s^2\tag{57}$$

where the first term on the right side of the equation reflects the viscous effects and the second term reflects the inertial effects. β is called the non-Darcy flow coefficient with a dimension of 1/length. This equation could be re-written in dimensionless form as:

$$f\_{pm} = \frac{-\frac{dP}{dx}}{\beta \rho V\_s^2} = 1 + \frac{1}{Re\_{pm}}\tag{58}$$

where *fpm* is the friction factor for flow in consolidated porous medium and *Repm* is the porous-medium Reynolds number, also referred to as Forchheimer number, defined as:

$$R c\_{pm} = \frac{\rho V\_s \beta k}{\mu} \tag{59}$$

At low Reynolds number (*Repm* < 0.1), the first term on the right side of the dimensionless Forchheimer equation, Equation (58), can be neglected and it reduces to Darcy's law. At high Reynolds number (*Repm* > 10), the second term on the right side of the dimensionless Forchheimer equation, Equation (58), can be neglected and it reduces to *fpm* = 1. Figure 8 shows the prediction of *fpm* as a function of *Repm* using the Forchheimer equation, Equation (58).

**Figure 8.** Prediction of porous medium friction factor as a function of porous medium Reynolds number. The data points are generated using the Forchheimer equation (Equation (58)).

In order to apply the Forchheimer equation, the value of non-Darcy flow coefficient is needed. It can be determined experimentally. The experimental procedure to determine β involves two steps: In the first step, the permeability *k* of the porous sample is determined from Darcy's law by restricting the experiments to low flow rates, and in the second step, the experiments are conducted at high flow rates and β is evaluated directly from the Forchheimer equation from the knowledge of pressure drop versus flow rate. The non-Darcy flow coefficient decreases with the increase in the permeability of the porous medium [20].

The Forchheimer equation is equally applicable to unconsolidated porous medium such as a packed bed of particles. Ergun equation is a special case of the Forchheimer equation. The Forchheimer equation, Equation (57) reduces to the Ergun equation, Equation (36), when permeability *k* and non-Darcy flow coefficient β are replaced by the following expressions:

$$k = \frac{d\_p^2 \varepsilon^3}{150(1 - \varepsilon)^2} \text{ and } \beta = \frac{7}{4} \left(\frac{1 - \varepsilon}{d\_p \varepsilon^3}\right) \tag{60}$$

The relationships between *fpm* and *fb*, and *Repm* and *Reb* are as follows:

$$f\_b = 1.75 \, f\_{pm} \tag{61}$$

$$\mathcal{Rc}\_b = \left(\frac{600}{7}\right) \mathcal{Rc}\_{\text{pm}} \tag{62}$$

#### 3.2.2. Entropy Generation in Flow through Consolidated Porous Media

Consider one-dimensional flow of an incompressible Newtonian fluid in consolidated porous medium, as shown in Figure 9. The rate of entropy generation in consolidated porous medium per unit volume is given by Equation (47), re-written as:

$$\dot{S}\_{\rm G}^{\prime\prime\prime} = \frac{V\_s}{T} \left( -\frac{dP}{d\mathbf{x}} \right) = \left( \frac{1}{T} \right) \left( \beta \rho \, V\_s^3 \right) f\_{\rm pm} \tag{63}$$

The porous medium friction factor *fpm* for laminar flow (*Repm* < 0.1) is given as:

$$f\_{pm} = \frac{1}{Re\_{pm}}\tag{64}$$

Consequently, Equation (63) reduces to:

$$
\dot{S}\_{G}^{\prime\prime\prime} = \left(\frac{1}{T}\right) \left(\frac{\mu}{k}\right) V\_{s}^{2} \tag{65}
$$

Thus, entropy generation rate per unit volume of the porous medium in steady laminar flow of a Newtonian fluid is directly proportional to the fluid viscosity, inversely proportional to the porous medium permeability, and directly proportional to the square of the superficial fluid velocity in the medium. The entropy generation rate also depends on the temperature.

**Figure 9.** One-dimensional flow in consolidated porous medium.

The porous medium friction factor in turbulent flow (*Repm* > 10) is given as:

$$f\_{pm} = 1\tag{66}$$

Consequently, Equation (63) reduces to:

$$
\dot{S}\_{G}^{\prime\prime\prime} = \left(\frac{1}{T}\right) (\rho \beta) V\_{s}^{3} \tag{67}
$$

Thus, entropy generation rate per unit volume of the porous medium in steady turbulent flow of a Newtonian fluid is independent of the fluid viscosity and is directly proportional to the cube of the superficial fluid velocity in the bed. The entropy generation rate also depends on the fluid density, the non-Darcy flow coefficient β, and the temperature.

To cover the full range of porous medium Reynolds number *Repm*, we substitute the Forchheimer equation, Equation (58), into Equation (63) to obtain the following equation valid over the full range of porous medium Reynolds number *Repm*:

$$\dot{S}\_{G}^{\prime\prime\prime} = \frac{1}{T} \left[ \left( \frac{\mu}{k} \right) V\_s^2 + (\rho \beta) V\_s^3 \right] \tag{68}$$

Equation (68) is simply the superposition of Equations (65) and (67).

#### **4. Simulation Results and Discussion**

#### *4.1. Entropy Generation in Flow through Packed Beds*

#### 4.1.1. Laminar Flow

The entropy generation rate in laminar flow through packed beds, per unit volume of the bed, as a function of superficial bed velocity is shown in Figure 10 for the following conditions: *T* = 298.15 K, *dp* <sup>=</sup> 1 mm, <sup>μ</sup> <sup>=</sup> 18.5 <sup>μ</sup>Pa·s, and <sup>ρ</sup> <sup>=</sup> 1.184 kg/m3. Equation (50) is used to generate the plots. For a given value of the bed porosity <sup>ε</sup>, the entropy generation rate . *S <sup>G</sup>* increases with the increase in the superficial velocity *Vs* of the fluid in the bed.

**Figure 10.** Effect of bed porosity on the rate of entropy generation per unit volume of the packed bed in laminar flow (*<sup>T</sup>* <sup>=</sup> 298.15 K, *dp* <sup>=</sup> 1 mm, <sup>μ</sup> <sup>=</sup> 18.5 <sup>μ</sup>Pa·s, <sup>ρ</sup> <sup>=</sup> 1.184 kg/m3). The data points are generated using Equation (50).

The increase in . *S <sup>G</sup>* with the increase in *Vs* is larger if the bed porosity ε is small. Additionally, the entropy generation rate falls sharply with the increase in the bed porosity. These results are as expected. With the increase in the fluid velocity, the rate of mechanical energy dissipation into frictional heat increases due to increases in viscous stresses and velocity gradients in the fluid. Consequently, there occurs an increase in the rate of entropy generation. When the bed porosity is increased, the bed structure becomes more open to fluid flow resulting in a decrease in the resistance to fluid motion and mechanical energy dissipation and, hence, a decrease in the rate of entropy generation.

Figure 11 shows the effect of fluid viscosity on the rate of entropy generation in laminar flow. The plots are generated using Equation (50). With the increase in fluid viscosity, the resistance to fluid motion increases due to an increase in the viscous stresses. With the increase in viscous stresses, the rate of mechanical energy dissipation into frictional heat (internal energy) increases. The conversion of highly ordered mechanical energy into disorderly internal energy is reflected in an increase in the rate of entropy generation.

**Figure 11.** Effect of fluid viscosity on the rate of entropy generation per unit volume of the packed bed in laminar flow (*T* = 298.15 K, *dp* = 1 mm, ε = 0.40, ρ = 1.184 kg/m3). The data points are generated using Equation (50).

The effect of particle diameter on the rate of entropy generation in laminar flow through packed beds is shown in Figure 12. A sharp reduction in the rate of entropy generation occurs with the increase in the particle diameter. With the increase in the particle diameter, the number density of particles is decreased for the same bed porosity ε. It can be readily shown that the number of particles per unit volume of bed *np* and the fluid-solids contact area per unit volume of bed *Asolids* are:

$$n\_p = \frac{\theta(1-\varepsilon)}{\pi d\_p^3};\ A\_{\text{solids}} = \pi d\_p^2 n\_p = \frac{\theta(1-\varepsilon)}{d\_p} \tag{69}$$

Equation (69) assumes that the particles are uniform spheres of diameter *dp*. The decrease in the number density of particles with the increase in particle diameter results in a decrease in contact area between the fluid and solids resulting in a decrease in the resistance to fluid motion. Consequently the rate of mechanical energy dissipation into internal energy and hence the rate of entropy generation decrease with the increase in particle diameter.

**Figure 12.** Effect of particle diameter on the rate of entropy generation per unit volume of the packed bed in laminar flow (*<sup>T</sup>* = 298.15 K, <sup>ε</sup> = 0.40, <sup>μ</sup> = 18.5 <sup>μ</sup>Pa·s, <sup>ρ</sup> = 1.184 kg/m3). The data points are generated using Equation (50).

Figure 13 shows the effect of temperature on the rate of entropy generation in laminar flow through packed beds. The rate of entropy generation decreases with the increase in temperature keeping other factors (fluid properties, bed characteristics) constant. At a given fluid velocity *Vs*, the rate of mechanical energy dissipation into internal energy is the same at different temperatures as the fluid properties are kept constant. Then why does the rate of entropy generation decrease with the increase in temperature? It so happens that the increase in entropy with the increase in internal energy is sensitive to temperature. For incompressible materials, the rate of change of entropy with respect to internal energy, that is, the derivative *dS*/*dU*, is given as:

$$\frac{dS}{d\varPi} = \frac{1}{T} \tag{70}$$

Higher the temperature, smaller is the increase in entropy for the same increase in internal energy.

**Figure 13.** Effect of temperature on the rate of entropy generation per unit volume of the packed bed in laminar flow (<sup>ε</sup> = 0.40, *dp* = 1 mm, <sup>μ</sup> = 18.5 <sup>μ</sup>Pa·s, <sup>ρ</sup> = 1.184 kg/m3). The data points are generated using Equation (50).

#### 4.1.2. Turbulent Flow

Figures 14–16 show the entropy generation rates for turbulent flow through packed beds under different conditions. Due to broad ranges of entropy generation rates and superficial fluid velocities, the plots are drawn using log-log scale. The plots on log-log scale are linear with a slope of 3 as expected from Equation (52). With the increases in bed porosity and particle diameter, the entropy generation rates decrease due to a decrease in the resistance to fluid motion and mechanical energy dissipation. There is no dependence of entropy generation on fluid viscosity in turbulent flow. However, fluid density now plays a role. With the increase in fluid density, the rate of entropy generation in turbulent flow increases due to an increase in the resistance to fluid motion and mechanical energy dissipation caused by inertial effects.

**Figure 14.** Effect of bed porosity on the rate of entropy generation per unit volume of the packed bed in turbulent flow (*<sup>T</sup>* = 298.15 K, *dp* = 1 mm, <sup>μ</sup> = 18.5 <sup>μ</sup>Pa·s, <sup>ρ</sup> = 1.184 kg/m3). The data points are generated using Equation (52).

**Figure 15.** Effect of particle diameter on the rate of entropy generation per unit volume of the packed bed in turbulent flow (*<sup>T</sup>* = 298.15 K, <sup>ε</sup> = 0.40, <sup>μ</sup> = 18.5 <sup>μ</sup>Pa·s, <sup>ρ</sup> = 1.184 kg/m3). The data points are generated using Equation (52).

**Figure 16.** Effect of fluid density on the rate of entropy generation per unit volume of the packed bed in turbulent flow (*T* = 298.15 K, ε = 0.40, *dp* = 1 mm, μ = 18.5 μPa·s). The data points are generated using Equation (52).

Figure 17 shows the entropy generation rate in a packed bed over a broad range of fluid superficial velocity covering the full range of packed bed Reynolds number. The plot is generated from Equation (53) under the following conditions: *T* = 298.15 K, ε = 0.40, *dp* = 1 mm, μ = 18.5 μPa·s, and ρ = 1.184 kg/m3. The predictions of Equation (53) overlap with the two asymptotes, Equation (50) for low velocities and Equation (52) for high velocities.

**Figure 17.** Prediction of . *S <sup>G</sup>* for packed beds over a broad range of fluid superficial velocity under the conditions: *<sup>T</sup>* <sup>=</sup> 298.15 K, <sup>ε</sup> <sup>=</sup> 0.40, *dp* <sup>=</sup> 1 mm, <sup>μ</sup> <sup>=</sup> 18.5 <sup>μ</sup>Pa·s, and <sup>ρ</sup> <sup>=</sup> 1.184 kg/m3. The data points are generated using Equation (53).

#### *4.2. Entropy Generation in Flow through Consolidated Porous Media*

The entropy generation rate in laminar flow through consolidated porous medium is given by Equation (65). The key factor affecting the . *S <sup>G</sup>* vs. *Vs* behavior is the ratio of fluid viscosity to permeability, <sup>μ</sup>/*k*. Figure <sup>18</sup> shows the plots of . *S <sup>G</sup>* vs. *Vs* for different values of μ/*k*. The plots on a log-log scale are linear with slopes of 2. With the increase in μ/*k* ratio, the entropy generation rate increases at any given superficial velocity *Vs*. When the fluid viscosity is increased (keeping *k* the same), the entropy generation rate increases due to an increase in mechanical energy dissipation caused by viscous stresses. When the permeability *k* of the porous medium is decreased (keeping μ the same), the porous medium becomes less permeable to fluid flow resulting in larger resistance to fluid motion and hence larger rates of mechanical energy dissipation and entropy generation. Note that *k* is typically in the range of 10−<sup>15</sup> to 10−<sup>12</sup> m<sup>2</sup> for porous sandstones [20].

**Figure 18.** Effect of viscosity to permeability ratio (μ/*k*) on the rate of entropy generation per unit volume of consolidated porous medium in laminar flow (*T* = 298.15 K, β = 10<sup>8</sup> m<sup>−</sup>1, ρ = 103 kg/m3). The data points are generated using Equation (65).

In turbulent flow through consolidated porous medium, the key factor governing the . *S <sup>G</sup>* vs. *Vs* behavior is ρβ where β is the non-Darcy flow coefficient. With the increase in ρβ, the entropy generation rate increases at any given *Vs*, as shown in Figure 19. Note that the non-Darcy flow coefficient is inversely related to the porous medium permeability [20]. For sandstones, β is typically in the range of 10<sup>8</sup> to 10<sup>12</sup> m−1. With the increase in β, the porous medium becomes less permeable to fluid flow resulting in an increase in the resistance to fluid motion. Consequently, the rates of mechanical energy dissipation and entropy generation increase with the increase in ρβ. The fluid density has a similar effect on entropy generation rate in the turbulent regime.

**Figure 19.** Effect of ρβ on the rate of entropy generation per unit volume of consolidated porous medium in turbulent regime (*<sup>T</sup>* <sup>=</sup> 298.15 K, <sup>μ</sup> <sup>=</sup> 10 mPa·s, *<sup>k</sup>* <sup>=</sup> <sup>10</sup>−<sup>12</sup> m2). The data points are generated using Equation (67).

Figure 20 shows the entropy generation rate in a consolidated porous medium over a broad range of fluid superficial velocity (10−<sup>4</sup> <sup>≤</sup> *Vs* <sup>≤</sup> 10 m/s) covering laminar, transition, and turbulent flow regimes. The plot is generated from Equation (68) under the following conditions: *T* = 298.15 K, <sup>μ</sup> <sup>=</sup> 1 mPa·s, <sup>ρ</sup> <sup>=</sup> 1000 kg/m3, *<sup>k</sup>* <sup>=</sup> <sup>10</sup>−<sup>12</sup> m2, <sup>β</sup> <sup>=</sup> 108 <sup>m</sup><sup>−</sup>1. The predictions of Equation (68) overlap with the limiting low *Repm* asymptote (Equation (65)) at low superficial velocities and with the high *Repm* asymptote (Equation (67)) at high superficial velocities.

**Figure 20.** Prediction of . *S <sup>G</sup>* for consolidated porous medium over a broad range of fluid superficial velocity under the conditions: *<sup>T</sup>* = 298.15 K, <sup>μ</sup> = 1mPa·s, <sup>ρ</sup> = 1000 kg/m3, *<sup>k</sup>* = 10−<sup>12</sup> m2, <sup>β</sup> = 10<sup>8</sup> m<sup>−</sup>1. The data points are generated using Equation (68).

#### **5. Conclusions**

In conclusion, a novel approach is described to teach the second law of thermodynamics via the analysis of flow through packed beds and consolidated porous media. The second law of thermodynamics and the relevant background in fluid mechanics are reviewed briefly. The link between entropy generation and mechanical energy dissipation in flow through packed beds and consolidated porous media is established in terms of the directly measureable pressure loss. Equations are developed to predict the entropy generation rates in a porous medium in terms of the flow variables, fluid properties, and structural properties of the medium. Simulation results dealing with entropy generation in a porous medium are presented and discussed. The proposed approach can be implemented at an undergraduate level either as an experimental exercise dealing with pressure loss measurements in flow through a porous medium such as a packed bed or as a simulation exercise. The material presented in the article is suited for third year engineering students after they have completed introductory courses in fluid mechanics and thermodynamics.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **CFD Julia: A Learning Module Structuring an Introductory Course on Computational Fluid Dynamics**

#### **Suraj Pawar and Omer San \***

School of Mechanical and Aerospace Engineering, Oklahoma State University, Stillwater, OK 74078, USA **\*** Correspondence: osan@okstate.edu; Tel.: +1-405-744-2457; Fax: +1-405-744-7873

Received: 1 July 2019; Accepted: 19 August 2019; Published: 23 August 2019

**Abstract:** CFD Julia is a programming module developed for senior undergraduate or graduate-level coursework which teaches the foundations of computational fluid dynamics (CFD). The module comprises several programs written in general-purpose programming language Julia designed for high-performance numerical analysis and computational science. The paper explains various concepts related to spatial and temporal discretization, explicit and implicit numerical schemes, multi-step numerical schemes, higher-order shock-capturing numerical methods, and iterative solvers in CFD. These concepts are illustrated using the linear convection equation, the inviscid Burgers equation, and the two-dimensional Poisson equation. The paper covers finite difference implementation for equations in both conservative and non-conservative form. The paper also includes the development of one-dimensional solver for Euler equations and demonstrate it for the Sod shock tube problem. We show the application of finite difference schemes for developing two-dimensional incompressible Navier-Stokes solvers with different boundary conditions applied to the lid-driven cavity and vortex-merger problems. At the end of this paper, we develop hybrid Arakawa-spectral solver and pseudo-spectral solver for two-dimensional incompressible Navier-Stokes equations. Additionally, we compare the computational performance of these minimalist fashion Navier-Stokes solvers written in Julia and Python.

**Keywords:** CFD; Julia; numerical analysis; finite difference; spectral methods; multigrid

#### **1. Introduction**

Coursework in computational fluid dynamics (CFD) usually starts with an introduction to discretization techniques of partial differential equations (PDEs), stability analysis of numerical methods, the order of convergence, iterative methods for elliptic equations, shock-capturing methods for compressible flows, and development of incompressible flow solver. A lot of emphases is put on the theory of discretization, stability analysis, and different formulations of governing equations for compressible and incompressible flows. Students can gain a better understanding of CFD subject with actual hands-on programming of numerical solutions to mathematical models that present different fluid behavior. To develop a flow solver from scratch is a fairly complex task. Therefore, almost all developmental CFD courses begin with one-dimensional problems to explain theory and fundamentals and then build on to complicated problems like compressible or incompressible flow solvers.

CFD Julia is a programming module that contains several codes for problems ranging from one-dimensional heat equation to two-dimensional Navier-Stokes incompressible flow solver. Some of these problems can be included as part of programming assignments or coursework projects. We include different types of problems, different techniques to solve the same problem and try to introduce CFD using a practical approach. There are a couple of programming modules available online related to CFD. CFD Python learning module is a set of Jupyter notebooks that tries to teach

CFD in twelve steps [1]. This module also contains bonus modules for numerical stability analysis, and advanced programming in Python using Numpy [2]. It starts with an introduction to Python and ends with the development of incompressible flow solver for lid-driven cavity problem and channel flow. There is an online module available for a course on numerical methods covering finite difference methods [3]. The material for this coursework was created using IPython notebooks. HyperPython is a set of IPython notebooks created to teach hyperbolic conservation laws [4]. HyperPython module covers high-resolution methods for compressible flows and also teaches how one can use the PyClaw package to solve compressible flow problems. PyClaw is a hyperbolic PDE solver in 1D, 2D, and 3D and it has several features like parallel implementation and choice of higher-order accurate numerical methods [5].

In this work, we use Julia programming language [6] as a tool to develop codes for fundamental CFD problems. Julia is a new programming language which first appeared in the year 2012 and it has gained popularity in the scientific community in recent years. There is a number of programming languages that are already available such as Fortran [7], Python [8], C/C++ [9], Matlab [10], etc. Many of the state of the art CFD codes are written in Fortran and C/C++ because of their high-performance characteristics. Python and Matlab programming languages have simple syntax and students usually use these languages with ease in their coursework. Since the developer's language is different from the user's language, there will always be a hindrance to developing numerical codes. Many times a user has to resort to developer's language like C/C++ for getting the fast performance or have to use other packages to exchange information between codes written in different languages. Julia tries to solve this two-language problem. Julia language was designed to have the syntax that is easy to understand and will also have fast computational performance. Julia shows that one can have machine performance without sacrificing human convenience. Julia tries to achieve the Fortran performance through the automatic translation of formulas into efficient executable code. It allows the user to write clear, high-level, generic code that resembles mathematical formulas. Julia's ability to combine high-performance with productivity makes it the perfect choice for scientists working in different scientific domains.

In this paper, we use finite difference discretization for different types of PDEs encountered in CFD. For all problems, we provide the main script in Julia so that reader will get familiar with the Julia syntax. We start with the heat equation as the prototype of parabolic PDE. We present an explicit forward in time numerical scheme and implicit Crank-Nicolson numerical scheme for the heat equation. We also cover the multistage Runge-Kutta numerical approach to show how we can get higher temporal accuracy by taking multiple steps between two time intervals. We derive an implicit compact Pade scheme for the heat equation. This will give the reader an understanding of high-resolution numerical methods which are widely used in the direct numerical simulation (DNS) and large eddy simulation (LES). All fundamental concepts related to finite difference numerical methods are covered in Section 2 with the help of heat equation. This will lay the foundation of CFD and will familiarize the reader with basic concepts of discretization methods which will be used further for more complex problems.

We then present the hyperbolic equation in Section 3 using the inviscid Burgers equation. Hyperbolic equations allow having a discontinuous solution. Therefore, higher-order numerical methods have been developed that can capture discontinuities such as shocks. We present weighted essentially non-oscillatory (WENO) formulation for finite difference schemes and show their capability to capture shock using one-dimensional case. We show the mathematical formulation of WENO schemes, their implementation in Julia for two different boundary conditions. We use Dirichlet boundary condition and periodic boundary condition for the inviscid Burgers equation. This will help the reader to understand how to treat different boundary conditions in CFD problems. We also present a compact reconstruction of WENO scheme which has better accuracy and resolution characteristic than the WENO scheme with the price of solving a tridigonal system. Furthermore, we present discretization of the inviscid Burgers equation in its conservative form using different finite difference

grid arrangement. This type of formulation is more applicable to finite volume method. We present a method of flux splitting and constructing flux at the interface using a Riemann solver approach in Section 4. We use a Riemann solver approach to solve the one-dimensional Euler equation and is presented in Section 5. The Euler equation solver is demonstrated for Sod shock tube problem using three different Riemann solvers: Rusanov scheme, Roe scheme, and Harten-Lax-van Leer Contact (HLLC) scheme.

We present different methods for solving elliptic partial differential equations using the Poisson equation as an example in Section 6. The Poisson equation is a boundary value problem and is encountered in solution to incompressible flow problems. We show different methods for solving the Poisson equation and their implementation in Julia. We demonstrate direct solver using fast Fourier transform (FFT) for periodic boundary condition problem and fast sine transform (FST) for Dirichlet boundary condition. We also present different iterative methods for solving elliptic equations. We demonstrate the implementation of Gauss-Seidel iterative method and the conjugate gradient method in Julia. The readers will learn the application of both stationary and non-stationary iterative methods. We also show the implementation of a V-cycle multigrid solver to accelerate the convergence of iterative methods.

We show solution to incompressible Navier-Stokes equation with vorticity-streamfunction formulation in Section 7. This is one of the simplest formulations to Navier-Stokes equation as there is no coupling between the velocity and pressure. This allows us to use collocated grid and not a staggered grid for the Navier-Stokes equation. We present a periodic boundary condition case using the vortex-merger problem. The Dirichlet boundary condition case is demonstrated using the lid-driven cavity problem. We use the direct solver to solve the elliptic equation giving the relation between vorticity and streamfunction. However, any of the iterative methods can also be used. We use second-order Arakawa numerical scheme which has better energy conservation properties for the discretization of the Jacobian term. The reader will get familiar with different types of numerical methods apart from standard finite difference schemes with these examples. This type of formulation can also be used for different types of problems such as two-dimensional decaying turbulence, and Taylor green vortex problem. Additionally, we present hybrid Arakawa-spectral solver for two-dimensional incompresible Navier-Stokes equations in Section 8. The flow solver is hybridized using explicit and implicit scheme for time integration. Also, the Navier-Stokes equations are solved in its Fourier space with the nonlinear term computed in physical space using Arakawa finite difference scheme and converted to Fourier space. Section 9 details the implementation fully pseudo-spectral solver for two-dimensional Navier-Stokes equations. We show two approaches to reduce aliasing error: 3/2 padding and 2/3 padding rules. We also show the comparison between CPU time for codes written in Julia and Python for solving the Navier-Stokes equations.

CFD Julia module covers several topics pertinent to both compressible and incompressible flows. This module provides different finite difference codes that will help students learn not just the fundamentals of CFD but also the implementation of numerical methods to solve CFD problems. Using these examples, students can go about developing their solvers for more challenging problems for their coursework or research. We make all these codes available on Github and are accessible to everyone (https://github.com/surajp92/CFD\_Julia). Table 1 outlines all codes available in the CFD Julia module with their Github index. Table 1 presents a summary of topics covered and numerical schemes employed for problems included in the CFD Julia module. We write results for all problems into a text file and then use separate plotting scripts to plot results. We provide installation instruction for Julia and required packages, and plotting scripts in Appendices A and B.



\* The codes for these problems are provided in Python in addition to Julia language.

#### **2. Heat Equation**

We start with the one-dimensional heat equation, which is one of the simplest and most basic models to introduce the concept of finite difference methods. The one-dimensional heat equation is given as

$$
\frac{\partial u}{\partial t} = a \frac{\partial^2 u}{\partial x^2},
\tag{1}
$$

where *u* is the field variable, *t* is the time variable, and *α* is the diffusivity of the medium. The heat equation describes the evolution of the field variable over time in a medium.

There are different methods to discretize the heat equation, such as, finite volume method [11], finite difference method [12], finite element method [13], spectral methods [14], etc. In this paper, we will study the finite difference method to discretize all partial differential equations (PDEs). The finite difference method is comparatively simpler and is easy to understand in the initial stage of learning CFD. The finite difference method converts the linear or nonlinear PDEs into a set of linear or nonlinear ordinary differential equations (ODEs). These ODEs can be solved using matrix algebra techniques.

Before starting with the derivation of finite difference scheme for the heat equation, we define a grid of points in the (*t*, *x*) plane. Let Δ*x* and Δ*t* be the grid spacing in space and time direction, respectively. Therefore, (*tn*, *xi*)=(*n*Δ*t*, *i*Δ*x*) for arbitrary integers *n* and *i* as shown in Figure 1. We use the notation *<sup>u</sup>*(*n*) *<sup>i</sup>* for *u*(*tn*, *xi*).

The finite difference method uses Taylor series expansion to derive the discrete approximation for PDEs. Taylor series gives the series expansion of a function *f* about a point. Taylor series expansion for a function *f* about the point *x*<sup>0</sup> is given by:

$$f(\mathbf{x}\_0 + \Delta \mathbf{x}) = f(\mathbf{x}\_0) + \frac{\Delta \mathbf{x}}{1!} \frac{\partial f(\mathbf{x}\_0)}{\partial \mathbf{x}} + \frac{\Delta \mathbf{x}^2}{2!} \frac{\partial^2 f(\mathbf{x}\_0)}{\partial \mathbf{x}^2} + \dots + \frac{\Delta \mathbf{x}^n}{n!} \frac{\partial^n f(\mathbf{x}\_0)}{\partial \mathbf{x}^n} + \mathcal{T}\_n(\mathbf{x}),\tag{2}$$

where *n*! denotes the factorial of *n*, and T*<sup>n</sup>* is a truncation term, giving the difference between the Taylor series expansion of degree *n* and the original function. The difference between the exact solution to the original differential equation and finite difference approximation is termed as the discretization or truncation error. The leading term in the discretization error determines the accuracy of the finite difference scheme. The finite difference scheme for any PDE must have two properties: consistency and stability to ensure the convergence. We would like to suggest a textbook "Finite difference schemes and partial differential equations" [15] which explains numerical behavior of finite difference scheme in detail. In this paper, we focus on the application of the finite difference method to study problems encountered in CFD.

**Figure 1.** Finite difference grid for one-dimensional problems.

#### *2.1. Forward Time Central Space (FTCS) Scheme*

We derive the forward time central space (FTCS) numerical scheme for approximating the heat equation. Using the Taylor series expansion in time for the discrete field *<sup>u</sup>*(*n*) *<sup>i</sup>* , we get

$$
\mu\_i^{(n+1)} = \mu\_i^{(n)} + \frac{\partial u\_i^{(n)} }{\partial t} \Delta t + \mathcal{O}(\Delta t^2), \tag{3}
$$

$$\frac{\partial u\_i^{(n)}}{\partial t} = \frac{u\_i^{(n+1)} - u\_i^{(n)}}{\Delta t} + \mathcal{O}(\Delta t). \tag{4}$$

Equation (4) is the first-order accurate expression for the time derivative term in the heat equation. To obtain the finite difference formula for second-derivative term in the heat equation, we write the Taylor series expansion for *<sup>u</sup>*(*n*) *<sup>i</sup>*+<sup>1</sup> and *<sup>u</sup>*(*n*) *<sup>i</sup>*−<sup>1</sup> about *<sup>u</sup>*(*n*) *i*

$$
\mu\_{i+1}^{(n)} = \mu\_i^{(n)} + \frac{\partial u\_i^{(n)}}{\partial x} \Delta x + \frac{\partial^2 u\_i^{(n)}}{\partial x^2} \frac{\Delta x^2}{2!} + \frac{\partial^3 u\_i^{(n)}}{\partial x^2} \frac{\Delta x^3}{3!} + \mathcal{O}(\Delta x^4), \tag{5}
$$

$$u\_{i-1}^{(n)} = u\_i^{(n)} - \frac{\partial u\_i^{(n)}}{\partial x} \Delta x + \frac{\partial^2 u\_i^{(n)}}{\partial x^2} \frac{\Delta x^2}{2!} - \frac{\partial^3 u\_i^{(n)}}{\partial x^2} \frac{\Delta x^3}{3!} + \mathcal{O}(\Delta x^4). \tag{6}$$

Adding Equations (5) and (6) and dividing by Δ*x*2, we get the second-order accurate formula for second-derivative term in heat equation

$$\frac{\partial^2 u\_i^{(n)} }{\partial \mathbf{x}^2} = \frac{u\_{i+1}^{(n)} - 2u\_i^{(n)} + u\_{i-1}^{(n)}}{\Delta \mathbf{x}^2} + \mathcal{O}(\Delta \mathbf{x}^2). \tag{7}$$

Substituting Equations (4) and (7) in Equation (1), we get the FTCS numerical scheme for the heat equation as given below

*Fluids* **2019**, *4*, 159

$$\frac{u\_i^{(n+1)} - u\_i^{(n)}}{\Delta t} = a \frac{u\_{i+1}^{(n)} - 2u\_i^{(n)} + u\_{i-1}^{(n)}}{\Delta x^2} \tag{8}$$

and we can re-write Equation (8) as an explicit update formula

$$
\mu\_i^{(n+1)} = \mu\_i^{(n)} + \frac{a\Delta t}{\Delta x^2} (\mu\_{i+1}^{(n)} - 2\mu\_i^{(n)} + \mu\_{i-1}^{(n)}).\tag{9}
$$

We denote the above formula as <sup>O</sup>(Δ*t*, <sup>Δ</sup>*x*2) formula meaning that the leading term of the discretization error is first-order accurate in time and second-order accurate in space. We can obtain the higher-order formula using the larger stencil to get better approximations. If the initial solution to the heat equation is given, then we can proceed in time to get the solution at future times using Equation (9). We can observe from Equation (9), that the solution at the next time step depends only on the solution at previous time steps. These numerical schemes are called an explicit numerical scheme. They are simple to code. However, explicit numerical schemes are restricted by some stability criteria. The stability criteria restricts the maximum step size we can use in spatial and temporal direction and hence explicit numerical schemes are computationally expensive.

We will now demonstrate how one can go about finding the solution at the final time step given an initial condition. We use the computational domain *<sup>x</sup>* <sup>∈</sup> [−1, 1] and *<sup>α</sup>* <sup>=</sup> 1/*π*2. The initial condition for our problem is *u*(*t* = 0, *x*) = −sin(*πx*). We use Δ*x* = 0.025 and Δ*t* = 0.0025 for spatial and temporal discretization. First, we initialize the array *u* and assign each element of the array using an initial condition. The solution after one time step is calculated using Equation (9) and we will store the solution in first column of matrix un[k,i] (variable in Julia) with the shape (*Nx* + <sup>1</sup>) × (*Nt* + <sup>1</sup>), where *Nx* is the number of divisions of the domain and *Nt* is the number of time steps between initial and final time. There is no need to store the solution at every time step, and we can just store the solution at intermediate time steps in a temporary array. Here, we store the solution at every time step to see how the field *u* evolves. The main script of the code which computes and stores the numerical solution at every time step is given in Listing 1.

We compare the solution at final time *t* = 1 using the analytical solution to the one-dimensional heat equation. The analytical solution for Equation (1) is given by

$$
\mu(t, \mathbf{x}) = -\mathcal{e}^{-t} \sin(\pi \mathbf{x}).\tag{10}
$$

We also compute the absolute error at *t* = 1 using the below formula

$$\epsilon(\mathbf{x}\_i) = |u\_i^{\text{exact}} - u\_i^{\text{numerical}}|. \tag{11}$$

**Listing 1.** Implementation of FTCS numerical scheme for heat equation in Julia. ✞ ☎

```