**Teaching and Learning of Fluid Mechanics, Volume II**

Editor

**Ashwin Vaidya**

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

*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: http://www.mdpi.com).

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**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1998-2 (Hbk) ISBN 978-3-0365-1999-9 (PDF)**

© 2021 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 Editor**

**Dr. Ashwin Vaidya** is a faculty member in the Department of Mathematics at Montclair State University in New Jersey, USA. He also serves as the director of the Complex Fluids Laboratory and holds a joint appointment in the Department of Physics and Astronomy at MSU. His research interests lie in the general areas of complex systems which includes fluid mechanics, fluid structure interaction, biological fluid mechanics and non-equilibrium thermodynamics, among others. He is also interested in issues surrounding math and science education with a focus on making 'creativity' a central theme of our educational practice.

**Ashwin Vaidya**

Department of Mathematics and Complex Fluids Laboratory, Montclair State University, Montclair, NJ 07043, USA; vaidyaa@mail.montclair.edu

#### **1. Introduction**

This issue showcases a compilation of papers on fluid mechanics (FM) education, covering different sub topics of the subject. The success of the first volume [1] prompted us to consider another follow-up special issue on the topic, which has also been very successful in garnering an impressive variety of submissions.

As a classical branch of science, the beauty and complexity of fluid dynamics cannot be overemphasized. This is an extremely well-studied subject which has now become a significant component of several major scientific disciplines ranging from aerospace engineering, astrophysics, atmospheric science (including climate modeling), biological and biomedical science and engineering, energy harvesting, oceanography, geophysical and environmental science and engineering, etc. While each of these disciplines has its own nuances and specific constraints, the fundamental physics behind the kinds of 'flow' phenomena discussed remains the same. In this volume, we bring together articles from authors with diverse expertise ranging from mathematics, physics, mechanical engineering, aerospace engineering, environmental engineering, and chemical engineering to discuss topics in fluid mechanics, many of which are of multidisciplinary interest.

The focus of all articles in this issue remains on the presentation of fundamental and advanced ideas on fluid mechanics which are suitable for presentation in an undergraduate or graduate course in fluid mechanics. Overall, I would divide the collection into the following four categories: (a) Pedagogy of fluid mechanics; (b) experimental or lab-based perspectives; (c) computational approaches; and (d) mathematical fluid mechanics. The following pages provide a brief summary of each of the contributions.

#### **2. Pedagogical Issues**

Student-centered practices such as problem-based and project-based learning (PBL) are more commonly practiced in the arts. PBL related instructional methods promote a more inductive approach to learning whereby generalizations and abstractions follow from first understanding specific cases. This approach is in contrast to the deductive strategy taken in the sciences which is a more top-down approach and a possible cause of alienation towards math and science in several students. The concept of problem-based learning began more than 30 years ago in the context of medical education and has been defined as the "posing of a complex problem to students to initiate the learning process" [2,3] and as "experiential learning organized around the investigation and resolution of messy, real-world problems." [4]. PBL can be implemented at various scales in a course with a focus from a "teacher to student-centered education with process-oriented methods of learning." [5]. The recent popularity of the project-based learning approach in physics and engineering education is based on research indicating the effectives of PBL in enhancing student engagement [5–7]. This volume presents a selection of papers that speak to the efficacy of PBL-based experiences in fluid mechanics courses.

The first article in the collection [8] by authors Garrard, Bangert, and Beck discusses an innovative pedagogical approach by planning and delivering large scale, multi-disciplinary labs with as many as 80 students in a single cohort and nearly 1000 students over a year.

**Citation:** Vaidya, A. Contributions to the Teaching and Learning of Fluid Mechanics. *Fluids* **2021**, *6*, 269. https://doi.org/10.3390/fluids6080269

Received: 23 July 2021 Accepted: 28 July 2021 Published: 30 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

Topics discussed in the lab are those which would be common to students from various branches of engineering at their institution including basic flow measurements, pressure driven flow in pipes, dimensional analysis, and other open ended projects. Besides the obvious financial and logistic advantages that come from sharing of resources, the authors point to its pedagogical value and efficiency.

The second article by Pérez-Sánchez and López-Jiménez [9] in this collection highlights the use of PBL in a fluid mechanics course taught at a Hydraulic Engineering Department in Spain that caters to over 2000 students. Just as in an earlier paper [8], the project described in this article proposes the coordination of fluid-based labs in different subjects at both the bachelor's and master's degree levels. The paper discusses the improvement in student performance, as well as the new teaching approaches which the authors note to have "increased the student's satisfaction index".

The final article in this section by Zoupidis, Spyrtou, Pnevmatikos, and Kariotoglou [10] is aimed towards teaching the concept of 'floating and sinking' (FS) to elementary school children. The authors explain the value of a "density-based explanatory model ... rather than the buoyancy-based" arguments typically used to explain FS phenomena, which is a conceptually challenging concept for children who instinctively associate floating and sinking with visceral experiences of 'lightness' and 'heaviness', respectively. The paper presents and evaluates the success of a novel instructional design paradigm founded on inquiry-based learning.

#### **3. Experimental in Fluid Mechanics**

The first article on experiments in FM by Wulandana discusses an impressive student driven project focused on building a recirculating flow tank [11]. Such tanks are an essential part of the collection of any fluid mechanics-based program and are very valuable due to their versatility and the ease with which many topics can be easily introduced. The only drawback of a prefabricated tunnel is its prohibitive cost. In this article, the author describes the design and fabrication of an open flow tank, built by students as part of their senior design project in a mechanical engineering program. The construction of this tank additionally provides opportunity for training in computer aided design (CAD) and computational fluid dynamics (CFD) as students perform comparative tests to validate flow structures in ideal and experimental conditions and also improve experimental designs so they meet expected flow conditions. The ideas introduced in this article can be replicated in any engineering program and also lead to interdisciplinary learning opportunities for students in physics and mathematics.

In the classic book, *A Splash of a Drop*, published by A.M. Worthington in 1895 [12], the world was introduced to the stunning visual world of droplet splashes. Worthington was certainly a strong influence on the movement in fluid dynamics scholarship to incorporate a visual element in order to understand the complex and beautiful structures that lie hidden behind the veil of transparency. Advancements and affordability of optical technologies makes is easier to introduce students to the fascinating world of flow visualization. The paper by Moghtadernejad, Lee, and Jadidi [13] introduces us to a course on multiphase flow where the instructors lead students on an experimental and theoretical investigation of splashes. Students also investigate the impact of temperature, wettability, impact velocity, droplet volume, shape, and relative humidity upon the splash dynamics. As the instructors note, this is an apt topic to introduce in a fluid mechanics course since it brings advanced knowledge into the classroom and also provides opportunities to discuss fundamental science and applications.

#### **4. Computational Approaches**

In this section we feature computation-based articles which cover both, articles of methodological nature and those that use computations to illustrate interesting physics of flows which can be introduced in any course on fluid mechanics.

The article [14] by Oz and Kara discuss computational methods relevant to 'Boundary Layer theory', a subject appropriate for introduction in an upper level undergraduate or graduate course. Relevant problems such as Blasius, Hiemenz, Homann, and Falkner–Skan flow equations are derived and numerically solved using the language, Julia. The codes have also been made freely available to the readers.

The article [15] by Ahmed, Pawar, and San provides a fundamental introduction to the mathematics and computational aspects of data assimilation methods which are fundamental to the study of climate science. Readers are exposed to the various methodologies through a series of Python modules which can be easily incorporated and adapted in an advanced course which treats such methods.

Mou, Wang, Wells, Xie, and Iliescu provide a survey of reduced order models which are computational models "whose dimension is significantly lower than those obtained through classical numerical discretizations" [16]. ROMs, in their various forms, have been found to be valuable in several complex computations involving uncertainty quantification, control, and shape optimization and in the numerical simulation of fluid flows. In this article, the authors summarize recent developments in ROM for barotropic vorticity equations, which are used to model geophysical flows.

In the article by Mongelli and Battista [17], the authors undertake a systematic study of pendulum dynamics by properly and fully accounting for the flow around a moving body which is not captured through the classical mechanical pendulum equations (see, also, another recent study that examines this issue through the lens of the least action principle [18]). The authors develop a computational fluid dynamics (CFD) model of a pendulum using the open-source fluid-structure interaction (FSI) software, IB2d. Comparisons with the results of the classical ODE model reveal very interesting and noteworthy results which ought to be discussed in any class which discusses pendulum dynamics.

Karlson, Nita, and Vaidya [19] discuss the interesting physics behind the vortex shedding phenomena. Computations using the program COMSOL are used to analyze the length of the primary vortex behind an elliptical body with varying eccentricities. The vortex length is used as a metric to understand and identify flow transitions from steady symmetric to asymmetric regimes which could potentially also be used as a noninvasive experimental strategy to distinguish flow regimes. The impact of the eccentricity of the body is seen to be particularly significant. While the physics itself is interesting and easy to follow and can even be discussed in an elementary FM course, such an example can be easily implemented in an advanced course in fluid mechanics or CFD course where students are exposed to a software for flow modeling.

#### **5. Mathematical Fluid Mechanics**

The final paper by Berselli and Spirito [20] is an extremely well written and much needed review of one of the most challenging mathematical problems of the last two centuries [21] and listed as one of the 'Millennium problems' in mathematics, namely the existence of solutions to the Navier–Stokes equation (NSE). While the history of this problem and various approaches is long and complex, the authors have done an excellent job in explaining and leading the readers through one aspect of this problem, namely the global existence of Leray–Hopf weak solutions to the NSE. I would strongly recommend that this article be made part of any course in an upper level undergraduate or even in an early graduate course in theoretical fluid mechanics or PDEs.

The papers in this volume, while selective and covering various different topics, showcase significant and cutting-edge knowledge of fluid mechanics in a manner that is easily adaptable for presentation in a course to undergrads, graduate students, and even in K12 settings. While the foundational materials traditionally taught in FM courses are important, our texts and curricula have not changed very much since the middle of the last century. The articles here provide templates for 'lesson plans' which can easily be implemented in our courses to make them more current and up-to-date. Many of the computation focused articles provide plug-and-play codes that can be implemented

without much training and time and comprising the larger objectives of the courses. We hope that educators will take note and find these papers helpful in their own teaching efforts and also in encouraging their own efforts towards incorporating other newer results into their classroom discussions.

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

#### **References**


## *Article* **Large-Scale, Multidisciplinary Laboratory Teaching of Fluid Mechanics**

#### **Andrew Garrard \*, Krys Bangert and Stephen Beck**

Multidisciplinary Engineering Education, University of Sheffield, Sheffield S10 2TN, UK; k.bangert@sheffield.ac.uk (K.B.); s.beck@sheffield.ac.uk (S.B.)

**\*** Correspondence: a.garrard@sheffield.ac.uk

Received: 18 October 2020; Accepted: 9 November 2020; Published: 11 November 2020

**Abstract:** The nature of fluid mechanics makes experimentation an important part of a course taught on the subject. Presented here is the application of a novel, large-scale multidisciplinary model of practical education in a fluids engineering laboratory. The advantages of this approach include efficiencies through the economy of scale leading to better pedagogy for students. The scale justifies dedicated academic resources to focus on developing laboratory classes and giving specific attention to designing activities that meet learning outcomes. Four examples of applying this approach to fluid mechanics experiments are discussed, illustrating tactics that have been developed and honed through many repeated instances of delivery. "The measurement lab" uses a flow measurement context to teach identifying and managing general experimental uncertainty. In this lab, new students, unfamiliar with fluid mechanics, are guided through a process to gain understanding that can be applied to all future experimental activities. The "pressure loss in pipes" lab discusses the advantage of and process for sharing equipment and teaching resources between multiple cohorts. Here, the provision for students is adapted for context, such as the degree program or year of study. The "weirs big and small" lab provides a methodology for teaching the power of dimensional analysis to mechanical engineers using a field of fluid mechanics that is outside their usual theoretical studies. Finally, the "spillway design" lab discusses mechanisms for delivering independent, open-ended student experiments at scale, without excessive staff resource requirements.

**Keywords:** practical engineering education; fluid mechanics; learning and teaching; laboratories

#### **1. Introduction**

Laboratory practicals are often included as part of the scheduled delivery for courses teaching physical principles. They allow students an opportunity to understand the physical manifestation of underlying concepts and compare theoretical models to real world results and can cater for alternative learning styles. These justifications are pertinent for courses in fluid mechanics. The nature of the subject often involves the understanding of qualitative or counterintuitive concepts that are best understood through a tactile experience. It can be argued that a visceral instinct for the behaviour of fluids can only be obtained with sufficient experience of its application in the real world. In addition, teaching fluid mechanics usually involves understanding concepts and models to predict the behaviour of a fluid for specific scenarios, for example, flow in a pipe or around a wing. As these scenarios become more complex, models rely increasingly on the use of empiricism in order to overcome the inability of analytical methods to model the flow. Indeed, compared to in other engineering disciplines, the requirement to introduce empirical correlations into predictive models occurs for even relatively simple physical systems, such a turbulent flow in a straight, horizontal pipe. Understanding the value of experimental affirmation and validation is critical for the development of well-rounded students studying fluid mechanics. This has been clearly voiced in a recent publication, where the

authors describe a blended approach to experimentation [1]. Within a university, each department is traditionally responsible for the delivery of laboratory teaching alongside a range of other teaching methods, such as lectures, tutorials, problem classes, design classes etc, which combine to provide students with courses in a particular subject discipline. The Faculty of Engineering at the University of Sheffield have adopted an alternative approach. One department, Multidisciplinary Engineering Education (MEE), is only responsible for the laboratory practicals of all 10 subject-specific degree programmes in the Faculty of Engineering, allowing the other departments increased time and resources to focus on classroom-based teaching methods and academic research. The volume of practical teaching delivered by MEE is consequently an order of magnitude greater than that which would typically be undertaken by departments offering individual degree programmes.

There are a number of advantages to this teaching model [2], such as the increased efficiency of infrastructure in common business processes, the reuse of teaching resources across similar activities and consistent experiences for students across the practical portion of their curriculum. The collective purchasing allows the possibility of buying many identical copies of experimental equipment and justifying the operation of large-capacity laboratories, with a related pedagogical advantage.

MEE is housed in the University of Sheffield's purpose-built Diamond building. Among the facilities is the Fluids Engineering Laboratory, which is used by Mechanical, Civil, Aerospace, Chemical, Bio and General Engineering students, in the order of 1000 students per year group. The laboratory has capacity for 80 students, typically working in groups of four, and is staffed by dedicated members of academic and technical staff. The laboratory is equipped with 20 identical copies of each of the pieces of experimental apparatus used for different aspects of fluid mechanics, including 20 hydraulic benches, on which a variety of internal flow experiments can be performed, and 20 wind tunnels, with which external flow experiments can be performed. The capacity provides three advantages of particular relevance to learning and teaching.


The Fluids Engineering Laboratory in the Diamond has been operating since 2015. The subsequent period of delivery has allowed a great deal of expertise to be developed in the teaching of practical fluid mechanics classes. Presented here are four examples of lab classes that have been honed through many repeated instances of running the activities.

#### **2. The Measurement Lab**

This laboratory is taught to all engineering students as one of the first tasks they perform when arriving at the University (typically in the first or second week). It is designed to equip students with a healthy scepticism for the results displayed on instrumentation and a toolkit for dealing with the uncertainty inherent in all forms of experimentation. A hydraulic bench is connected to a fluidic circuit containing a series of flow measurement devices. Although, at this stage in their programmes, students will be unfamiliar with almost all fluid mechanics concepts, including those that underpin flow measurement, it is explicitly stated and reinforced throughout the teaching that the activity is about understanding and managing general experimental uncertainty. Part of the intention is to imply the universality of error and uncertainty of techniques for any experimental set up, even ones for which the concepts under investigation or the outcomes are unknown.

Despite the lack of technical understanding of fluid mechanics principles, it is reasonably straightforward to explain the concepts of the conservation of mass and, for an incompressible flow, conservation of the volume flow rate. Students are aware that they will be studying fluid mechanics as part of the engineering programmes and are keen to understand these basic concepts early as well as be introduced to real-world instrumentation that they may not have previously been exposed to in schools. As water passes from one device to another in series through the hydraulic circuit, it is evident to the students that the volume flow rate through each device must be identical. The students are tasked with predicting if the various devices will all record identical readings for the flow rate.

Prior to starting the activity, students are given a briefing and watch an instructional video discussing the methods for capturing uncertainty for various pieces of instrumentation, how to record it and how it can be propagated when raw data are processed. The activity involves applying these principles to unfamiliar equipment.

Students record raw data from the instrumentation: the heights from water columns attached to a Venturi meter and an orifice plate, the flow rate from a calibrated rotameter measuring in litres per minute and the timing of water collection using a measuring tank. In order to compare the measured flow rates, the raw data from each device need to be converted. This requires students to consider that, in order to compare, the same measuring unit is required for each device and provides an opportunity to discuss the relative merits of the more commonly used litres/minute over the S.I. unit of meters cubed per second. While UK students are often trained to habitually convert into S.I. units, there is little justification for doing so in this case, particularly as it eliminates the need to convert the results of the rotameter.

Students need to process the height difference between the water columns of the Venturi meter and the orifice plate into a flow rate, the physics of which will not yet be known to them. For the Venturi meter, a pre-prepared spreadsheet is provided, which outputs the flow rate when the water column heights are input. The spreadsheet solves Equation (1), which is provided to students to allow understanding of the mathematical relationship, but the process of manual calculation is not required. The spreadsheet provides the opportunity to repeatedly calculate different answers quickly and determine the relative impact of the uncertainty for different parameters. The square root function within this equation makes this process more interesting.

$$Q = C\_d A\_1 \sqrt{\frac{2g\Delta h}{\left(\frac{A\_1}{A\_2}\right)^2 - 1}}\tag{1}$$

A calibration chart, shown in Figure 1, is provided to convert the raw data from the orifice plate. The process of reading and extrapolating results from the graph introduces additional uncertainty in processing the raw data, and students are encouraged to consider how to best incorporate this into their calculated flow rates. This facilitates an opportunity for students to become comfortable with a limited availability of precision in a calculation, which is a common occurrence in the application of practical engineering.

**Figure 1.** Calibration chart for an orifice plate provided to students.

Students are guided through a process of calculating flow rates and the associated uncertainty of their results for the four flow measurement devices and then through a process of presenting this information graphically with an introduction to the concept of error bars. If conducted correctly, the results show that all the devices will record different values of the flow rate but, when the error bars are considered, all the results overlap within a certain region. Further discussion can be had about the methods for reducing uncertainty in the raw data and the advantages of in-line flow measurement compared to the volume displacement of the measuring tank.

This lab has been designed to achieve the specific learning outcomes of introducing the concepts of, methods to record and process for handling error and uncertainty in experimentation. Students are clearly made aware of this expectation and achieve these explicit learning outcomes as a result of participating in the activity. The same learning outcomes could have been achieved with a paper-based exercise, delivered outside a laboratory. However, the act of learning through doing is more likely to result in the concepts being retained by students and provides a real-world context in which to apply these skills in an engaging form that enhances the student experience.

As this activity is delivered to all undergraduate engineers, typically, in excess of 1000 students per year, the investment of time to develop high-quality instructional material and training teaching assistants is easily justified and makes the activity very resource efficient.

#### **3. Pressure Loss in Pipes Lab**

The pressure loss in horizontal pipes is measured in an experiment run in the Fluids Engineering Lab for Mechanical, Aerospace, Civil and Chemical Engineering students. This is an important part of the engineering curriculum [6]. The reuse of teaching material and equipment for multiple cohorts results in efficient resource utilization. However, teaching material is adapted and contextualized for specific degree programmes. Subject-specific nomenclature or units should be used appropriately for different engineering disciplines. For example, civil engineers would measure, record and process pressure in the units of meters of head, whereas aerospace engineers would typically use Pascals. In addition, the level of academic rigour and expectation for the students is adapted depending on the placement of the activity within the degree programme. At the University of Sheffield, mechanical engineering students perform the experiment in the first year, and civil engineering students, in the second. For civil engineers, less prescriptive instructions are provided, and more independence is expected while conducting the experiment.

During the experiment, the students will collect raw data from manometers to measure the pressure drop and a measuring tank to measure the flow rate. A range of different-diameter and roughness pipes are available on each of the 20 hydraulic benches. The raw data are processed into Reynolds numbers and empirically derived friction factors, allowing students to generate their own Moody diagram that can be compared to a published version. The objective of the laboratory experiment is not to develop expertise in performing the mathematics. Students are provided with a "guided calculation", where the steps to process the raw data into a processed result are described in the instructions and executed by students on one piece of data, to ensure they understand the mathematical methods. Breaking each part of the calculation into defined steps makes the debugging of errors, by the students or teaching assistants, more straightforward, which is necessary when dealing with large class sizes. Once students have demonstrated they understand the process, a spreadsheet to automate the calculations on the remaining data points is released.

Mandating students to complete the hand calculation before using the spreadsheet opens the opportunity for discussion with students that perceive the activity to be a trivial task only required to access the spreadsheet. Performing the hand calculation and using a tool to perform the calculations allows a two-way validation of each process, by comparing the results from each. When using any tool that has been provided, it is wise to ensure it operates as expected. Articulating the general merit of a validation approach, and how it can be applied to a student's future engineering tasks, can be used to place value on performing the task. In addition, digital collection allows an individual student's data to be pooled into a larger dataset that can be shared with the cohort, for the purposes of error and reproducibility analysis.

Having multiple benches to support a large class size presents an opportunity to improve efficiency beyond the economy-of-scale issues previously described. In this experiment, to determine the influence of the pipe specimen (diameter and roughness) on the friction factor, multiple pipes should be investigated. With one hydraulic bench or a small number of hydraulic benches, pipes need to be installed and removed to test the full range. With as many benches as specimens, benches can be set up with particular specimens, and students can move around the laboratory to each piece of apparatus. Learning the procedure, executing it and performing the subsequent bleeding of air all consume student time and cognitive capacity in ways that do not directly relate to the intended learning.

As the students come into a session well prepared and the lectures on the subject are fresh in their minds, they do not require much assistance to conduct the laboratory. Hence, it can be taught by four staff (Academic, Technical and Teaching Assistants). These staff are able to spend the time discussing the work and providing feedback to the students, resulting in a much richer experience for everyone.

#### **4. Weirs Big and Small Lab**

One of the hardest fluid mechanics topics for students to understand is the importance and power of dimensional analysis. This is because there are a number of difficult concepts when contemplating scales, which can be simplified through the correct application of dimensionless numbers. These include the fact that experiments are needed to be able to obtain the constants for every geometry in a given situation. An understanding of dimensional analysis is needed for studying aerodynamics and heat transfer, but the teaching of it generally suffers from two deficiencies: firstly, students typically become very tangled up in the details of the subject rather than the method application, and secondly, almost all of the work and examples involve Reynolds numbers (as in the Pressure Loss in Pipes lab described above). Thus, an interesting approach to illuminating this topic has to have two requirements: not using Reynolds number and extracting some constants that are then applied to different scales of equipment. If it could be engaging, challenging and fun as well, that would be even better.

The hydraulic benches can be configured to allow students to perform open channel flow experiments with sharp-edged weir plates (square and triangular) and can measure the water height over the weir and water flow rate. The cohort for whom the lab was created are second year Mechanical Engineers whose curriculum does not contain free surface flow. This presents an ideal opportunity to introduce this topic to these students while showing the power of dimensionless groups. Students are shown, using the tools from their lectures (Buckingham π theory), that the dimensionless groups involved in this type of flow are the Froude number √*<sup>V</sup> gH* and length ratio *<sup>H</sup> <sup>b</sup>* , as shown in Figure 2,

and that the volumetric flow rate (*Q*, m3·s<sup>−</sup>1) can be derived from Equation (2):

$$Q = \mathbb{C}\_D b \sqrt{\mathbb{g}} H^{\frac{3}{2}} \tag{2}$$

where *CD* is the discharge coefficient and must be experimentally ascertained for the given geometry. The objective of this experiment is to experimentally determine the discharge coefficient for the small weirs and see how this scales to a geometrically equivalent, larger weir that is installed in the lab's 10 m flume. The main learning outcome of this activity is for students to be able to see both the power of dimensionless numbers and how extracting the constants experimentally is a required part of the process.

**Figure 2.** Dimensions used in weir calculations and the experimental apparatus.

In order to optimise the time students spend in the laboratory, a comprehensive pre-experimental activity was created, material from which can be seen in Figure 3. This consisted of a series of presentations. These were either recorded with overheads and voiceovers for theory (the flow over weirs, which Mechanical Engineers do not cover as part of their course). Another online lecture on fitting exponentials to a series of *x* and *y* data was created to help students understand one way of turning experimental data into equations, which they need to do to extract the constants in this experiment. Short quizzes created in the Virtual Learning Environment using adaptive release ensure that the students engage with each presentation prior to moving on to a subsequent section. The students are presented with three videos on using the large flume, operating the small flow rigs and reading the Vernier scale on the large flume. This culminates in a quiz on reading the Vernier and finally a compulsory test on Health and Safety issues relevant to the experiment and laboratory space. The use of adaptive release means that students cannot get to the final test without completing all the previous ones.

**Figure 3.** Examples of instructional material, including pre-experimental videos, provided to students.

This preparation works well. Students move on very quickly to the practical experiment. As required, they calculate the *CD* on the small rigs from the flow and height measurements they recorded. For each session, the laboratory leader sets a different flow rate for the large flume, and the students predict the height above the weir that the water should reach before they measure it. This, in effect, gamifies the lab, as they are ascertaining their own experimental accuracy and are able to compare it with their peers'. They are thus able to grasp the power and value of dimensionless numbers and scaling in engineering. It shows them that methods such as the Buckingham π theory are merely tools to be used to solve real problems. It also illuminates the way that researchers need to use experiments to be able to create a set of results to identify and extract generalities. It is also a good way

to use our facilities, as they only perform a single experiment on the large, unique rig but can perform a large number on the multiple small rigs.

By the end of the experiment, it was noted that the students appreciated how it was possible to scale between different sizes of models, but how one needed to be aware of all the important parameters. In a more general sense, the intellectual understanding that this journey provides through models and sizes acts as a conceptual bridge towards a better understanding of the dimensionless numbers that they will use in their fluid mechanics and thermodynamics careers and education.

#### **5. Spillway Design Lab**

Teaching large cohorts presents a tension: investing significant resources in activities that can be reused by many students is efficient and provides a high-quality, professional laboratory experience, but prescribed activities can limit the opportunity for students to explore open-ended activities and, for example, learn through failure. This tension can be partially overcome with the application of a multidisciplinary approach. Within MEE's portfolio of practical engineering education is manufacturing and fabrication. This provides a holistic integration of making that is available to staff and students.

With significant manufacturing capability, department-based workshop staff and tools, in-house builds of bespoke teaching equipment are feasible. Typically, engineering teaching equipment for use with students would have been bought from suppliers. The two significant downsides of this approach are that it is extremely expensive compared to an in-house build (if full-time staff time for design, fabrication and prototyping is excluded) and the equipment is not designed to achieve specific learning outcomes.

Unlike their Mechanical counterparts, second year Civil Engineering students study open channel flow and the design of flow control devices during their second semester. As they are reasonably advanced students, having been prescriptively taught the fundamentals of operating in a laboratory environment in preceding years, their practical activities are designed to be conducted independently, open ended and genuinely experimental, i.e., conducting empirical work to discover something previously unknown. To achieve these outcomes, a bespoke experimental rig was conceived, as shown in Figure 4.

**Figure 4.** Bespoke rig for teaching spillway design.

The teaching of design principles is outside the remit of a fluid mechanics course, but a holistic approach to a programme can be used to apply previously learned design principles to the design of a weir and spillway to fulfil specific problem specifications. As such, the bespoke rig was designed to provide a constrained number of parameters that can be adjusted. Designated adjustable parts can be fabricated by students using readily available equipment, such as laser cutters or 3D printers, and inserted into the rig for testing in the fluids lab flume. Students are expected to design their adjustable parts based on theory delivered during lectures, predict the flow and test their predictions in the laboratory.

Prior to being given access to the flume, students are provided with extensive equipment and Health and Safety training to be allowed access to the laboratory without staff supervision. Compliance is established with online tests, and keys to the room/equipment are provided by reception staff who check for completion of the test. Students are able to book use of the equipment at a time convenient to them, and academic staff time input is minimised. This approach provides students with an opportunity to exert agency over their own learning.

Without a dedicated team of staff focused on providing students with an exemplary practical experience and the multidisciplinary team of academics and technicians working collaboratively, the development of practical teaching and bespoke equipment that is unobtainable from suppliers would be significantly more difficult for departments to justify resourcing.

#### **6. Discussion and Conclusions**

Historically, fluid mechanics laboratories have been run as a type of cottage industry, with lecturers specifying and delivering a couple of labs across the year on a single piece or possibly a couple of pieces of equipment. This meant that students could receive this laboratory at any time over a year. Each experiment had to be free standing; many students would conduct the experiment long before or after they were introduced to the theory, missing a crucial window for learning reinforcement. This meant that different students would, in effect, obtain different learning outcomes depending on their understanding of the background to the topic. The staff members responsible for the delivery of the lectures, tutorials, exams etc. set their own labs; they tended to be similar in difficulty, scope and assessment (usually a report). There was no coherence or progression along the course and, without the capacity to focus exclusively on the laboratory activities, very little in the way of designing teaching with constructive alignment towards the overall learning objectives.

MEE's multidisciplinary approach of professionalising and integrating the practical experience of the students allows many of these common issues to be obviated. The scale of the laboratories allows entire cohorts to perform an experiment in a short time period so that practical and theoretical work can be interwoven and used to support each other. In many cases, the formal lecture becomes the introduction to the laboratory. Many Electrical Engineering departments have rooms set up with multiples of equipment, but this approach is rare outside engineering. Having dedicated staff who deliver the only practical experience to a cohort, it is possible to *curate* an entire, progressive student experience starting from the closed and didactic (such as the Measurement Lab) and progressing to open-ended investigations such as the Spillway Design Lab. The result is that students receive an integrated and progressive learning experience culminating, after their first two years, in them becoming capable, reflective and autonomous experimenters ready to start independent project work.

As well as the efficient use of space and staff time, the experiments form a portfolio of work that can be renewed and repurposed as and when required. For example, within weeks of the creation (and delivery) of the Weirs laboratory, a lecturer from Civil Engineering asked if there was a laboratory for open channel flow for their MSc students. Not only was the answer "yes", but a version of all of the teaching and introductory material was ready for use. The laboratory sheet only needed updating to reflect the different approach to the theory and nomenclature used by a different discipline, but this was a minor investment of time and allowed the students to have an excellent practical experience to support their learning that would have been impossible under a different organisation.

There are, however, two potential drawbacks to the multidisciplinary approach, but these can be ameliorated if properly anticipated. Firstly, when the practical and theoretical teaching on a single module is delivered independently by different members of staff in different departments, the experience and messaging received by the students could become disconnected and incoherent. MEE overcomes these issues by setting up communication channels between the academics delivering classroom and practical teaching, allowing them to agree on how the labs are presented within the context of a module and ensure that the messaging to students is consistent. The tactics for achieving this include using material presented in lectures as part of laboratory tuition and vice versa. Secondly, there is a requirement for strong leadership within the faculty. The multidisciplinary model will only work if all the departments benefiting from the service agree to contribute to its resourcing. With any shared resource, issues of perceived value and equity for contributors can cause tension if not carefully managed.

Thus, in conclusion, there are a number of major advantages to teaching at scale in fluids laboratories, such as the efficiency, temporal proximity to lectures and scalability. Due to the integration and professionalisation of the practical teaching, it allows an integrated, progressive approach to student practical skills development to be implemented. Progressing from the usual method of teaching practical fluid mechanics to the new one demonstrated in the examples above is a difficult, long and potentially extremely expensive journey. We hope that we have shown you that the outcomes from it are worthwhile.

**Author Contributions:** Conceptualization, S.B., A.G. and K.B.; Methodology, S.B., A.G. and K.B.; writing—review and editing, S.B., A.G. and K.B. All authors have read and agreed to the published version of the manuscript.

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

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

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 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/).

## *Article* **Continuous Project-Based Learning in Fluid Mechanics and Hydraulic Engineering Subjects for Di**ff**erent Degrees**

#### **Modesto Pérez-Sánchez \* and P. Amparo López-Jiménez**

Hydraulic and Environmental Engineering Department, Universitat Politècnica de València, 46022 Valencia, Spain; palopez@upv.es

**\*** Correspondence: mopesan1@upv.es; Tel.: +34-96-387700 (ext. 28440)

Received: 6 May 2020; Accepted: 13 June 2020; Published: 15 June 2020

**Abstract:** Subjects related to fluid mechanics for hydraulic engineers ought to be delivered in interesting and active modes. New methods should be introduced to improve the learning students' abilities in the different courses of the Bachelor's and Master's degree. Related to active learning methods, a continuous project-based learning experience is described in this research. This manuscript shows the developed learning methodology, which was included on different levels at Universitat Politècnica de València. The main research goal is to show the active learning methods used to evaluate both skills competences (e.g., "Design and Project") and specific competences of the students. The research shows a particular developed innovation teaching project, which was developed by lecturers and professors of the Hydraulic Engineering Department, since 2016. This project proposed coordination in different subjects that were taught in different courses of the Bachelor's and Master's degrees, in which 2200 students participated. This coordination improved the acquisition of the learning results, as well as the new teaching methods increased the student's satisfaction index.

**Keywords:** outcomes competences; hydraulic engineering; hydraulic teaching; active methodology

#### **1. Evolution of Teaching in the University**

#### *1.1. New Paradigms in Hydraulic Engineering Teaching*

Hydraulics disciplines are in a higher number of degrees related to engineering topics [1]. Civil and environmental engineering courses are an example, although they are not exclusive. There are different Bachelors' and Masters' degrees, in which the fluid mechanics and hydraulic topics are present inside of the students' curricula, such as the mandatory or optative subject.

Teaching involves many methods to reach the learning results. Some of them are: master courses (i.e., a theoretical lesson taught by a professor); design projects; practical activities in the hydraulic lab; and informatic sessions, among others. All actions must provide students an integral and continuous vision of the hydraulic engineering (from fluid mechanics to environmental problems). However, the new students must experiment a necessary change in the new learning methodologies, allowing the students to reach the professional competences satisfactorily. Currently, the European Higher Education plans to decrease the credits to teach, increasing the required skills acquisition by activities, which are not an on-site class [2,3].

#### *1.2. The Significance of the Learning Habilities*

One of the problems that students must face is the complexity of numerous concepts in hydraulics subjects (e.g., fluid mechanics). The lecturer usually teaches theoretical matters and the student has to reach the learning results (e.g., master course, lectures, and exercises) using the professors' information (e.g., bibliography and exercises). Therefore, students must learn materials on their own with minimum guidance by professors. This methodology had good results in the last decades [4]. However, newer student generation demands the development of new teaching methods that rely on new tools. The use of active learning methods is highly recommendable since the students participate in the learning process actively [5]. These methods are based on student activities (e.g., 'playing and learning' using simulations, project-based learning, and role activities) they are a possible solution to improve the students' learning. Continuous project-based learning was proposed across different levels in 2016 as a part of these active strategies [6,7]. The current research shows the results, which were obtained by the coordination between bachelor and master matters from 2016 to 2018.

Currently, numerous researches show the professors should introduce the new learning tools and activities (e.g., simulations, experimental cases, and playing learning) using information and communication technology [8]. Using these tools engages students actively in the learning process. In this line, the Universitat Politècnica de Valencia (UPV) carries out the 'UPV generic students' outcomes 2015–2020 [9]. The main goal of this project is the introduction of 13 generic outcomes, which will improve the students' skills and their curricula [9].

To adapt the new strategic plans, an innovation and educational improvement project has been implemented between different professors of the Hydraulic and Environmental Engineering Department of the UPV since 2016. The main objective is to establish a transversal and vertical coordination in different subjects. The purpose is the acquisition improvement of the learning results by the students. Therefore, the project proposes an evaluation methodology using different rubrics according to the domain level (depending on the year and degree). Besides, the research compares the different subjects in different courses.

In this particular case, this proposal allowed students to start a hydraulic project draft in fluid mechanics topics (e.g., students sized a water branched network). Furthermore, they continued its development in hydraulic machinery matters (e.g., students designed a pump system) and finally, they designed a total project in their last matter fluid facilities (e.g., students sized fluid facilities in a hotel). The development of hydraulic projects, as a learning strategy, is a methodology that was proposed in other universities some years ago. When this method is used, the students must plan, implement, and evaluate complete works, which are applied in real case studies [8,10].

Project-based learning (PBL) allows students to acquire key knowledge and skills through the development of projects that respond to real-life problems [11]. The objective is to enhance students' autonomy. They become the main actor of their own learning process. This training evolves introducing new complex tasks each course using the same project [8]. In this learning process, professors guide and support the students throughout the entire project.

#### *1.3. Hydraulic Engineering Learning Challenges*

The stage of the studies of hydraulic engineering must be in constant evolution since the future professionals must face great challenges. These are aligned on terms of sustainability and optimization of the management. Therefore, hydraulics subjects at the university level cover many fields such as: urban hydraulics, watershed management, the pollutants dispersion, hydraulic machinery, river dynamics and restoration, water resources management, hydraulic works, aspects of flows to sheet free of charge involved in sanitation, the water-energy nexus and many other subjects that are being taught in different faculties masterfully. In all cases, the hydraulic engineering is present in the core subject, such as fluids mechanic and/or hydraulic machinery in the engineering bachelor's degree (e.g., electrical, mechanical, and chemistry).

Currently, the knowledge transfer requires the future students must be autonomous and capable. They have to develop skills, which allow them to solve the new challenges. The present manuscript shows the continuous project-based learning experience, which has been developing at UPV. This practice is focused on hydraulics subjects, which are teaching both the bachelors' and masters' degree. The proposed teaching project increases the development of real projects, decreasing the hours of lessons. This time decrease is complemented with online material (e.g., teaching video and laboratory tutorials) including workshops and specific conferences. These sessions are developed by companies or guest speakers in the university focused on students.

This research is a good example for engineering bachelors' and masters' degrees related to hydraulic and environmental topics at different levels. The manuscript summarizes teaching methodology and results, which developed a teaching project. The experience was carried out at Universitat Politècnica de València. Two thousand and two hundred students participated in thirteen hydraulic subjects, which were part of the teaching project and they were from different years. The students worked the hydraulic concepts using a methodology, in which they reached the learning results through the development of hydraulic projects. The strategy enabled to evaluate both specific and outcomes competences. Before this teaching project, the students were not evaluated of their skill competences and they did not use active methods. Previous to this project, the students' training was based on master courses and laboratory practices. The participation was up to 80% and the student's satisfaction was measured by surveys.

#### **2. Materials and Methods**

#### *2.1. Structure of the Hydraulic Engineering for a Student of a Bachelor's and Master's Degree in the UPV*

When the structure of the hydraulic engineering was analyzed at the UPV, there was a complete interweaving with other matters in their different bachelor's and master's degrees. These subjects (Figure 1), which were distributed throughout student training, were: (i) basis subjects in hydraulics and fluid mechanics; (ii) subjects related to hydraulic machines; (iii) subjects related to hydroelectric plants and wind power machinery; (iv) subjects related to hydraulic facilities; (v) materials in oleo hydraulic and pneumatic systems; (vi) matters in relation to the water-energy binomial; (vii) matters in computational fluid dynamics (CFD) modeling; (viii) matters in relation to the hydraulic aspects of wastewater treatment; and (ix) matters in relation to the dispersion of contaminants in receiving fluid media.

**Figure 1.** Structure hydraulic engineering subjects at Universitat Politècnica de Valencia (UPV; BD.—Bachelor's degree and MD.—Master's degree).

#### *2.2. Learning Proposal Based on Learning Projects at Di*ff*erent Levels, Developing the Transversal Competence "Desing and Project"*

The project-based learning (PBL) is a methodology focused on learning, research, and reflection. In this methodology, the students should reach the correct solution of a problem, using an autonomous and continuous learning. This problem was proposed by the lecturer once he/she teaches the theoretical concepts [12]. This methodology was included on the teaching project in the hydraulic and environmental engineering department [13]. It was applied on different matters, which were part of different courses and levels (i.e., Bachelor's and Master's). Therefore, the teaching project got a continuous project-based learning (CPBL) in the students' training [14]. Besides, the CPBL application at different training times of the student enables one to work and evaluate different transversal competences (e.g., time planning, permanent self-learning, and oral communication as well as design and project (DP)).

Time planning was proposed for each subject and it must be followed by both students and professors through the different phases. These steps (Figure 2) were divided on face-to-face and non-face-to-face lessons. The first phase allows students to know the theoretical concepts throughout the master course, the development of computer practices as well as the development of basic problems related to the taught issue. Once these are known, the learning results of each unit should be practiced in progressive development of the project. This practice is non-face-to-face and the students must use information from the UPV webpage. In this section, they have supplementary material. Along this phase, the students work on self-learning and the professor gives them help in group meetings. The collaboration between the professor and students improves the acquisition of the learning results.

**Figure 2.** Example of temporal distribution to reach the learning results in a subject (U is unit, T is a test, and SP is simulation practice).

The different activities, their planning, and their dedication were defined using good practices sheets. This sheet was developed by lecturers, and it enabled to coordinate the subject between students and professors, and coordination improved when different teachers participated in teaching the subject.

#### *2.3. Proposal of Rubrics*

One of the objectives to develop an active learning is the definition of evaluation items. In this particular case, different rubrics were introduced to combine the evaluation of the student's competences (i.e., skills and concepts). These rubrics were composed of different indicators, which had

four different descriptors for each one. The indicators measured the acquisition degree of the learning results. Table 1 shows the proposed indicators, which were used to do the proposed rubric (Appendix A). This rubric was used to evaluate the project in hydraulic machines. Other rubrics were published for different subjects: the wastewater network project [6], fluid mechanics [7], or fluid facilities in the chemistry industry [15]. The used rubrics were different in Master's and Bachelor's degrees, considering the reached level in the descriptor. In the Bachelor's degree, the pupils had to design a project with a level of draft. In contrast, the indicators were higher in the Master's degree, as the developed project should be more specific, and students must be more autonomous.


**Table 1.** Definition of weighted in the different indicators.

The attached rubric in Appendix A shows the new proposed rubric for hydraulic machines in which specific and skills competences were evaluated. The symbiosis between specific and transversal competences was developed using a matrix, which contained weights and ponderations. The first discrimination was done between 'not done' and 'developed task but the minimum is not reached'. When the student did not develop the descriptor, the numeric value was zero. If the student did the task, but it was not reached, the considered value was 3. If the descriptor was C, the numeric value was 5. When the descriptor was B, the considered value was 7 while the value was 10 when the descriptor was A. Each indicator had a specific weighted value, which was justified in Table 1.

Therefore, the numeric mark (NM) of the specific competence was obtained using Equation (1):

$$\text{NM} = 0.05I\_1 + 0.075I\_2 + 0.20I\_3 + 0.50\frac{\left(\sum\_{l=1}^{1-7} I\_4\right)}{7} + 0.075I\_5 + 0.10I\_6 \tag{1}$$

This expression enabled one to get the numeric value through descriptors. If the NM was less than 4.5, the learning result was "No reach—D" for transversal competence. When the value was between 4.5 and 6, the learning result was "In Development—C". A learning result of "Good—B" was reaching, when the NM was between 6 and 8. Finally, 'Excellent—A' was reached when the NM was greater than 8.

#### **3. Results**

#### *3.1. Students, Subjects, and Proposal of Projects*

The project was developed between 2016 and 2019, although it continues currently. In these years, thirteen subjects were taught in Bachelor's (second, third, and fourth year) and Master's (industrial engineering and hydraulic and environmental engineering) degree at UPV. One thousand and one hundred students participated and 13 professors from the hydraulic and environmental engineering department collaborated in the project each year. The manuscript shows results for all subjects although main subjects of the students' curricula (i.e., fluid mechanics, hydraulic machines and fluid facilities) were described deeply in this research. Fluid mechanics was taught in a second-year course of the Bachelor's Degree. This subject was in chemical, electrical, and mechanical engineering degrees (in this particular case, the results were related to the mechanic engineering degree). The students were between 19 and 25 years. One hundred and seventy students were involved, who were divided into two groups for theoretical teaching and four groups for practical classes.

Hydraulic machines was taught in the third-year course of the mechanical engineering degree. The subject was focused on analyzing the pumps-operation principles (velocities triangle and Euler's equation) as well as the machines selection and their regulation according to demand. The students were between 20 and 26 years. One hundred and forty students were involved. These pupils were divided into two theory groups and four practical groups.

Fluid facilities was taught in the first-year course of the industrial engineering Master's degree, after students achieved their Bachelor's degree. The subject contained the analysis of the different types of the fluid facilities, that is, water distribution networks, gas networks, and waste-water networks as well as the facilities, which involved the comfort in society (ventilation and hot sanitary water). For each facility type, the normative, design, analysis, and regulation were analyzed, applying it to the real cases study. The students were between 23 and 30 years. Three hundred and fifty students were involved. These students were divided into seven theory groups and twenty-one groups in practical lessons.

In relation to the fluid mechanics subject, an elemental water supply network was proposed, in which students proposed different diameters for pipes, considering flow and pressure conditions. Once the system was sized, the students had to analyze it using Epanet software [16], considering the constrain conditions (e.g., demand, level node, minimum pressure, and maximum velocity). The sizing was developed as a function on demand over time, using the uniform hydraulic slope criterion. The students run an extended period simulation, analyzing the pressure and flow variations in the different lines. Finally, the students proposed a short budget, considering both length and the chosen material. In this case, the project was supervised by the lecturer. The initial information, which was available for students was: network topology, reservoir head, water demand in each point over time, modulation curve for the different consumption patterns, and the minimum operational conditions of the network as well as the material type and cost of the pipelines. The work was focused on establishing a methodology to develop hydraulic calculus, encompassing the Bernoulli's and continuity equations. The students compared the different studied scenarios as a function of demand pattern using Epanet software.

The evaluation was a formative type. The students did meetings with the professor and they show partial results. The professor verified the calculus and solutions, proposing improvements to students. There were two meetings. The first meeting included the proposal of the network. The second meeting addressed the sizing of the water system. The correction of the project was developed using rubric (Appendix A), and a third meeting was done to explain to students the errors in the project.

The proposed work in hydraulic machines was individual. The activity was focused on analyzing the energy consumption and regulation of a pumped system. This water network was supplied considering two options. Option A: the water network was supplied from a reservoir, which was filled using a pump station and Option B: the water was directly supplied using pump systems. Option A enabled one to analyze the influence of the reservoir volume in the pump selection (mainly pumped flow) when the energy cost was considered (i.e., schedule and operation time). Option B was focused on applying the similarity laws, regulating the operation curve. The students had to define the rotational speed of the machine as a function of the demanded flow over time. The students defined: the control rules, the operation costs, and the efficiency parameters for each pump system, trying to minimize the cost per cubic meter. The student only had two constrains: demand over time and energy cost, which was the current Spanish energy price.

Finally, when the student undertook the fluid facilities subject, the proposed work had a higher level than previous tasks developed in the Bachelor's degree. At this time, the students were more mature and the cases were near real buildings. Therefore, their training should be more intensive and closer to reality. In 2017, the proposed activity was to develop a complete project (summarize, calculus supplements, drawings and budget, defining the qualities, and normative for the different used materials). The project was related to a complex building (e.g., hotel, hospital, and school since for each students' group it is different). In this project, the students had to connect basic knowledge of fluid mechanics and hydraulic machinery with the new learning results, which are reached in the fluid facilities subject. The students designed the different pipelines and equipment, which were necessary to supply the building (e.g., cold and sanitary hot water system, pumps, and ventilation, among others). This work was developed by teams, composed of three or four students. Once the work was finished, the students had to explain it in an oral session.

In all cases, the students' doubts were attended by teachers. Generally, the questions were solved by face-to-face meetings. However, the doubts solution was also solved using mail and/or a video conference. Throughout the process, the student contacted the professor to validate the different items of the project in each one of the phases and stages.

#### *3.2. Analysis of Results*

#### 3.2.1. Results

Figure 3a shows marks distribution in a students' group for hydraulic machines. Each indicator value can be observed for each student. The project mark was the upper 8/10 for 24 students while there were only six students who qualified below 5/10. Each indicator (from I1 to I6) is described in Table 2 and they are drawn in the Figure 3a.

Figure 3a shows the students worked really well I3 and I4 indicators. These were focused on the development of the simulations and the establishment of the control rules in the pumped systems to guarantee the hydraulic constrains (i.e., flow and pressure). In contrast, I5 was the worst developed indicator and it focused on the analysis and discussion of the results. However, the results were highly satisfactory.

Figure 3b,c shows the transposition from the mark to transversal competence in the different subjects that participated in the teaching project (Table 2). If observing the topic hydraulic machines (12659), 77% of students reached the A and B descriptors when the "Design and Project" competence was evaluated. Similar results were obtained in the rest of subjects shown in Figure 3b.

If all subjects were observed the satisfaction was higher, considering all students who participated in the teaching project. The participation in the project development was 82%, considering there were 1051 students in thirteen different matters in 2018 (1149 students were in 2017). When the "Design and project" competence was evaluated, 361 students reached an excellent degree (A). The B degree was reached by 286 students while 154 and 58 (6.75%) students obtained a C and D degree, respectively (Figure 3b).

Figure 4 shows there is a lineal relationship between exam and project marks in two years (i.e., 2016/2017 and 2017/2018). Therefore, the development of the activity helped students to acquire the hydraulic concepts as well as the methodology. Although there were no exams when PBL was applied, in order to compare the previous (traditional method using master courses) and new methodology (CPBL), an exam was proposed. This improvement contributed to reaching the learning results favorably. This trend was observed in majority of the studied subjects. Besides, when the project delivery was after the exam, the test mark did not have a relationship between them. Therefore, there was a greater significance to establish the date delivery before the test. The final marks were compared with previous years. The score increased around the 1–2 point about 10, reducing the number of students who failed the subject (6% in 2016/2017 and 8% in 2017/2018).

**Figure 3.** (**a**) Indicator in hydraulic machines related to Table 2; (**b**) results of the descriptors in the "Design and Project" competence for the subjects of the Bachelor's degree defined in Table 3; and (**c**) results of the descriptors in the "Design and Project" competence for the subjects of the Master's degree defined in Table 3.

Figure 4 shows the correlation between project and exam marks was strongly correlated when the students did not do the project correctly or they got a mark up to six. However, when the students developed an excellent project (mark between 7 and 9), they did not always get an excellent mark in their exams. It could be due to the student's collaboration and working on other skills such as analysis and resolution of problems. When they did the exam, the help between partners was not there, and therefore, they had to solve their doubts, which in some cases were not resolved correctly.





#### 3.2.2. Surveys

Two different surveys were developed for each subject. First survey gave information related to the subject planning. This survey helped to analyze if the coordination between taught concepts and project development was correct. Related to this, five questions were proposed (Table 3). These questions were related to: (i) the application of the activity with the concepts, which were taught in the classroom (Q1); (ii) the synchronism between activity and taught concepts (Q2); (iii) if the index developed by the professor to explain the methodology was clear (Q3); (iv) if the project development helps student to reach the learning results and train to do the evaluations (Q4); and (v) if the student would be interested in development of a project considering different subjects of the degree (Q5).

**Figure 4.** Relationship between the exam and project marks in hydraulic machines.

Figure 5 shows the results in the survey when it was done in hydraulic machines. The figure shows the results once 103 students (76%) answered it. There were 90% of students that positively agreed with Q1. This percentage was higher compared with other subjects in the UPV. This was a goal of the project, since it wanted to develop activities to increase the satisfaction in the students. These activities were focused on: (i) increasing the simulation lessons with software, (ii) visiting some buildings where the students can identify the studied facilities, and (iii) increasing the number of online videos in which they can visualize real solved case studies. In both years, the answer was similar between students. Therefore, they considered positive the use of this methodology to apply the teaching concepts.

**Figure 5.** Survey to analyze the development of the project in the topic: hydraulic machines for the 2017 and 2018 years.

The rest of the questions were mostly approved. They showed the majority of students agreed to develop the teaching methodology. Similar results were obtained in fluid facilities (Figure 6). The surveys analysis showed the student accepted this methodology although they had to invest more effort in the subject continuously. This learning obligated students to develop a planning to reach the objectives. In this case, results from only one year were presented, since the Master's students only undertook a course in the active methodology in 2017 (previously, they undertook a course in for their Bachelor's degree on the topic hydraulic machines using this teaching project).

**Figure 6.** Survey to analyze the development of the project in fluid facilities.

Figure 7 shows the results of a survey, which had four questions. The survey was proposed to students in the second-year or third-year level (once the student studied fluid mechanics). The questions measured the vertical coordination between subjects. The questions (Table 4) were related to: (i) the developed project that helped students to improve the acquisition of competences (Q6); (ii) if the development of the project, which was developed on fluid mechanics in the previous year, helped to improve the development of the project in hydraulic machines (Q7); if the previous study of the hydraulic concepts helped students to develop the project (Q8); if the use of a similar methodology between the project developed both fluid mechanics and hydraulic machines that helped students to develop the project in the hydraulic machines topic (Q9).

**Figure 7.** Survey to analyze the vertical coordination between subjects that are located on different courses and levels (Bachelor's or Master's degree).


**Table 4.** Questions related to coordination between years.

If Figure 7 was analyzed, it shows that the majority of students (upper 60%) considered the application of this methodology positively and it had influence on achieving good results in the development to their competences.

The developed experience verified that the students improved the acquisition of the learning results in the different subjects when they were compared with the previous years. Therefore, the professors' experience joined to the students' opinion show the development of active methodologies increased the positive attitude of the students. This emotional state made the students show a greater interest in the subject, improving their efficiency. However, this effect cannot occur in some cases, in which students think they learn a lot in these scenarios, but when tested they really are not. This occurs when the students do not work in the activities continuously and correctly throughout the year.

Similar strategies are being developed currently in UPV and other universities to motivate students to develop continuous learning. Currently, the development of projects is being planned at an institutional level. The learning project includes subjects that are part of different years and are in different areas. This situation improves the integration of the subjects in the students' curricula. Besides, the students understand better the significance of the different subjects when there is a global learning project. It occurs even though the matters are studied in different years. The project existence allows students not to view the subjects individually, interrelating the different matters.

This methodology can be extrapolated to other knowledge areas or degrees, adapting the projects to the learning results of each subject. The success of this methodology is verified in other countries and universities [2,3,17]. The development of the good practices sheet [16] and the definition of the learning goals allow one to organize the active methodology for any subject.

#### **4. Conclusions**

A case study was described in this research, which joined different hydraulic engineering topics. The subjects were taught using continuous project-based learning. The implementation of this methodology was new at UPV to develop the skills competences in the students. The development of the methodology from basic subjects (i.e., fluid mechanics) enabled one to define the procedure, which can be applied on subjects at the upper level. The methodology allowed professors to establish a schedule in which face-to-face time and a non-face-to-face lesson fit perfectly. Therefore, the use of a good practice sheet allowed students and professors to know their activities for each time. The use of these sheets improved the synchronization of the teaching (i.e., theoretical concepts, practices lessons, and activities) between them. Besides, the development of the good practice sheet helped professors to organize subject learning. The good practice sheet contained different tasks, which should be carried out by professors and students, defining the data and time of their development.

The proposed methodology is crucial to give students an action strategy when they have to develop similar projects in matters of the hydraulic area. The strategy improves the vertical coordination in the Bachelor's or Master's degree. This organization maximized the reach of the learning results since it mixed a face-to-face class and online videos and material as well as real projects to apply the taught concepts according to the students' capacity.

The methodology included a rubric for each subject. These evaluation criteria enabled us to evaluate the acquisition of the learning results, joining both specific and transversal competence of the 'design and project'. The rubric, which was used in hydraulic machines was shown in this manuscript (Appendix A). It defined the specific indicators and the descriptors, which are necessary to develop the project. Besides, the used expression, which correlated the specific and transversal competences in the students' curricula, was presented.

Two surveys were proposed to students. These questions showed the students' satisfaction for the structure of the activity. Besides, they considered it necessary to improve the acquisition of the learning results. The students were grateful of the use of this methodology in other subjects related to hydraulic engineering topics.

Finally, the new challenge in teaching should be focused on:


**Author Contributions:** The author P.A.L.-J. wrote the introduction and she analyzed the results of both tests and rubrics. The author M.P.-S. proposed the methodology and he contributed with the analysis of the results. All authors have read and agreed to the published version of the manuscript.

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

**Acknowledgments:** The described experience has been carried out as part of the work in the Innovation and Quality Education Teaching (EICE DESMAHIA) "Development of active methodologies and evaluation strategies applied to the field of Hydraulic Engineering" in the Universitat Politècnica de València.

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

#### **Appendix A**


**Table A1.** Specific indicators used to evaluate hydraulic machines.


#### **Table A1.** *Cont.*

#### **References**


© 2020 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/).

## *Essay* **Teaching and Learning Floating and Sinking: Didactic Transformation in a Density-Based Approach**

**Anastasios Zoupidis 1,\*, Anna Spyrtou 2, Dimitrios Pnevmatikos <sup>2</sup> and Petros Kariotoglou <sup>3</sup>**


**Abstract:** This essay synthesizes more than a decade of research, most of which has been published, on the teaching and learning of floating and sinking (FS) phenomena. The research is comprised of the iterative design, development, implementation and evaluation of a Teaching-Learning sequence (TLS) for the teaching and learning of density within FS phenomena. It was initiated within the frame of the European Community supported "Materials Science" project. Due to the many, different aspects of the project, each publication has focused on a particular part of the study (e.g., effectiveness and the iteration process). The didactic transformation for the teaching of FS phenomena is presented and discussed here. In doing so, it is essential to mention: (a) the students' ideas as the main cause of the scientific knowledge transformation, (b) the scientific/reference knowledge, and (c) the knowledge to be taught and its limitations. Thus, we intend to describe and justify the didactic transformation process and briefly synthesize the published (from previous papers) and unpublished results to show its effectiveness.

**Keywords:** inquiry-based instruction; science education; teaching-learning sequences; didactic transformation; primary level

#### **1. Introduction**

School children are familiar with floating and sinking (FS) [1], which is a main topic in the teaching of fluids in science education [2,3], especially at the primary and lowersecondary levels (10- to 15-year-olds). Although the topic is very common, and children have many everyday life experiences in FS phenomena, their interpretation is challenging, not only because of the difficulty of the scientific concepts and the respective explanatory models that are involved (e.g., density, buoyancy), but also because of these everyday experiences that students have and their subsequent ideas [4].

Research on FS has been extensive in the last few decades, both regarding students' ideas [5] and, consequently, about ways to effectively teach this topic [6]. Concurrently, Teaching-learning sequences (TLSs), i.e., medium-level curriculum unit packages, that include well-researched teaching-learning activities empirically adapted to student reasoning [7,8], are increasingly present in science education research, because they provide the opportunity to integrate teaching and learning theories and approaches, students' ideas about science concepts and explanations of natural phenomena, as well as the historical development of scientific concepts [7–10].

One of the most critical issues in the design and development of a TLS is the didactic transformation of the content, i.e., transforming the scientific knowledge into appropriate knowledge to be taught [7,11,12]. The choice of content and how it is transformed in order for it to be easily understood and readily adopted by students is crucial in the entire process of TLS development. Although this often takes place, none or very little of it is

**Citation:** Zoupidis, A.; Spyrtou, A.; Pnevmatikos, D.; Kariotoglou, P. Teaching and Learning Floating and Sinking: Didactic Transformation in a Density-Based Approach. *Fluids* **2021**, *6*, 158. https://doi.org/10.3390/ fluids6040158

Academic Editor: Ashwin Vaidya

Received: 26 March 2021 Accepted: 13 April 2021 Published: 14 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

usually conveyed. In other words, even though the didactic transformation process is an essential aspect in every teaching effort, and especially in the design and development of a TLS, researchers rarely describe the process in an explicit and detailed manner, possibly due to space restrictions [8,13]. Thus, colleagues who wish to further investigate any such didactic proposal's effectiveness do not have all the necessary information to repeat its implementation.

To describe the didactic transformation process of certain content, in our case, FS phenomena, it is important to mention, among other things: (a) students' ideas about the phenomenon and the concepts related to it, as the main reason for the scientific knowledge transformation, (b) the scientific/reference knowledge, and (c) the knowledge to be taught, in its new form, following scientific knowledge transformation, including its limitations.

Students' alternative ideas have played a decisive role in the planning of teaching in science education in the last few decades [6]. Consequently, students' ideas about FS phenomena, and the difficulties they face in adapting interpretations to be consistent with the scientific ones, need to be taken into account in every teaching effort that is developed within the frame of the prevalent constructivist approach [14]. Moreover, the study of the historical development of scientific knowledge concerning FS interpretations could contribute to the didactic transformation process by revealing the difficulties scientists had come up against in understanding and interpreting those phenomena throughout the centuries [14–16].

The TLS entitled "Density of materials in floating and sinking phenomena: Experimental procedures and modelling", which was initiated in the framework of the European Community supported "Materials Science" project (FP6, SAS6-CT-2006-042942), has been described elsewhere [17–22]. However, because there were many different aspects of content to be taught, i.e., declarative (density and floating sinking), procedural (control of variables strategy), and epistemological knowledge (nature and role of models), the focus of the previous published papers has been other than describing and justifying the didactic transformation, which we hope to do here.

In this paper, we briefly describe elements of our developmental research that have already been published, focusing, however, on the didactic transformation of content, as this has not yet been thoroughly presented or discussed and which we consider to be of paramount importance. Therefore, our aim is to describe the didactic transformation of the content of the TLS concerning FS phenomena, to underline the factors which influenced the process of its development, and to present the limitations of the transformation. Specifically, we justify the reasons why we chose the density-based explanatory model for FS phenomena, rather than the buoyancy-based (see Section 3), as well as providing arguments for the didactic transformation of the concept of density, which is still an open issue for an effective approach to FS learning. Furthermore, selected essential aspects of the revised version of the TLS, which was adapted from an initial study, are also described [19]. Moreover, a short presentation of both our published (from previous papers) [20,21], and unpublished results in FS [22], on the implementation of the revised TLS in a real-class environment is given. In this sense, we consider that this work is an original sample of a developmental research description in the framework of Design-Based Research approach in science education [8,23,24].

#### **2. Alternative Ideas and Difficulties in Explaining FS Phenomena**

Students seem to perceive FS phenomena visually. That is to say that they decide whether an object is in a floating or sinking state based on the object's position relative to the surface of the liquid [1]. For instance, the majority of students in Joung's study [1] answered that an object was floating in the water when at least a part of it was above the surface, most of whom chose the case where the object touches the surface of the liquid and fewer chose the case where the object was half-submerged. Also, the majority of students considered that an object had sunk in the water, in the cases where it was below the surface, i.e., (a) at the bottom, or (b) in the middle (between the bottom and the surface of the liquid), with a decreasing frequency of occurrence, respectively. In the same research, students considered an object just below the liquid surface either to be floating or to have sunk. Few students recognized that when the object was between the surface of the liquid and the bottom, then it is in a state between floating and sinking, and therefore, was suspended and remaining at rest in the liquid at the same location [1,25].

In addition, it seems that students explain and describe the phenomena in relation to perception-based macroscopic natural properties, such as weight, length, and volume [5,26–28]. In other words, students formulate their estimation concerning the floating of solid objects in a liquid by taking into account: (a) the heaviness/size of the objects, (b) the existence of hollows, (c) the existence of holes, (d) the interface/edge, orientation, shape and/or texture of the floating object, (e) the dimensions of the tanks in which floating takes place, (f) the amount and/or depth of the liquid, and (g) the liquid stickiness [29]. Needless to say, one of the most prevalent alternative ideas that students hold is that of case (a), that is, students most often claim that an object floats because it is small and/or light, and it sinks because it is big and/or heavy [26,30].

Consequently, when interpreting FS phenomena, students tend to focus on the properties of the objects or the liquids. Additionally, they seem to merely use causal linear reasoning, referring only to an object's or a liquid's property, instead of causal relational reasoning, which involves comparing object and liquid densities in their interpretations [31]. However, this is not the only obstacle in students using the specific causal relational reasoning to explain FS phenomena. Researchers who have studied students' conceptions of density [26,32,33] have found that they had difficulty in understanding this abstract concept. Firstly, students find it hard to understand the ratio of two quantities [34], such as that of mass per volume, particularly when those quantities are changing simultaneously [35]. Secondly, the concept of density is a property that is not directly perceived through the senses but can only be understood through mental reasoning and/or calculations [33,36]. Thirdly, students' difficulty in understanding density is rooted precisely in an already developed conceptual framework about matter and material kind [37], which is composed of perception-based physical quantities where the raw scientific notions of weight, volume and density coexist undifferentiated [33]. Consequently, these students consider density to be proportional to the size of an object or the object's quantity of matter.

To fully understand the reasons why an object floats or sinks, one needs to comprehend that the concept of water pressure is an intensive property, while the concept of buoyancy is a force, and not, as is usually the case with students, a property of an object, within the framework of Newtonian mechanics [4,5]. However, students very often confuse the FS states with the explanatory model; that is, they equate buoyancy, a construct/force in an explanatory model, with the state of floating. They also seem to think that buoyancy is a property of an object opposed to the interaction between an object and its surrounding fluid, as they are unable to understand buoyancy as a force, i.e., the interaction between two entities. Furthermore, many students have the misconception that the buoyancy of an object is inversely proportional to its density, while others are not sure about the direction of buoyancy [4].

In sum, students confront severe difficulties in interpreting FS phenomena in the framework of both density-based and buoyancy-based explanatory models. The reasons for this difficulty, however, are not the same in both cases. In the former, the difficulty is mainly due to the non-differentiation of the concepts of weight, mass, and density, in contrast to the latter case, where it is mainly due to the students' inability to understand the concept of force as the interaction between two entities.

#### **3. Floating and Sinking Teaching Approaches**

The way educators approach the teaching of FS phenomena can be put into two broad categories, according to the central concept of the explanation of the phenomena (Figure 1). In the first category are those cases that provide density-based explanations, e.g., [29], following the so-called elimination of variables approach. This approach focuses

on highlighting the variables that affect the FS phenomena in order to derive a prediction rule that will determine "which" body will float. In the second category are those cases that provide buoyancy-based explanations, e.g., [38], following the so-called scientific approach, that is, an interpretation using an equilibrium mechanism in order to explain "how" an object floats [15]. Several researchers provide both explanations concurrently [4,32].

**Figure 1.** Categories of teaching approaches to floating and sinking phenomena; recreated by the authors according to [15].

The scientific approach is rooted in Archimedes' and Galileo's explanatory models [15,39]. In brief, Archimedes explained the floating of a solid in an infinite container, using only the concept of weight as a quality/property of objects [40] and comparing the weight of the solid to the weight of the fluid that was displaced by the immersed portion of the solid. He proposed that the weight of the fluid displaced would be equal to the weight of the solid. However, Archimedes' model was effective only for FS phenomena in infinite containers [40]. On the other hand, Galileo overcame Archimedes' inadequacy by explaining floating of a solid in a finite container. Galileo's model comprised important discriminations of new-defined concepts [15,40], such as the distinction of floating and surface tension phenomena [41] and the differentiation of the concepts of weight and specific weight in a qualitative manner [15]. Both models derive their floatation laws from an analysis of the conditions that result in equilibrium. However, Galileo's model consists of more delicate and abstract concepts and, subsequently, can explain a larger range of phenomena [15,39]. On the other hand, the elimination of variables approach is rooted in Inhelder and Piaget's work [42]. These researchers were the first to record children's explanations on the floating phenomenon, focusing mainly on the ability of children to (a) classify a set of objects according to whether they float or sink in the water, and (b) explain the criterion by which the classification was made [15]. Moreover, they were interested in testing children's ability to eliminate inconsistencies in their initial explanations, such as using the weight of objects in order to interpret FS phenomena, and at the same time, to formulate the predictive floating law; that is, objects float if their density is less than that of water, where density is defined as the ratio of weight to volume. Thus, contrary to the scientific approach, which focuses on the construction of an interpretive model that will explain "how" a body floats, the elimination of variables approach focuses on the construction of a predictive model for FS phenomena, with the latter being less complicated, and subsequently, an easier process for students [15].

Usually, the density-based approach is adopted by primary schools and junior high schools [6,15], while senior high schools and colleges/universities adopt the buoyancybased approach [15] or a combination of both [4]. When one develops a curriculum for this topic, the chosen approach, whether density- or buoyancy-based, would also involve selecting different steps, practices and concepts in the teaching and learning process. Therefore, in the case of the elimination of variables approach being chosen, there would need to be an intermediate goal of the learning process, which would involve differentiating weight, volume, and density as this approach requires a clear understanding of these the concepts (Figure 1). Another intermediate goal would be explaining FS phenomena by using a partially correct explanatory model, that is, that an object's density as a property of materials determines its sinking or floating.

A correct explanatory model within the frame of the elimination of variables approach would be to use the concept of density in causal relational reasoning, that is, comparing the densities of an object and a liquid in order to come to a decision about the FS situation [31]. Contrastingly, in those curricula where FS is taught through an analysis of the equilibrium approach, it is not needed for students to discover or to be introduced to any intermediate/precursor concepts [29]. For example, sinking and floating can be explained as a result of the balance between gravity and buoyancy. In that case, the forces of buoyancy and gravity are the required concepts for the final learning goal. However, both buoyancy and gravity are scientific—rather than intermediate/precursor—concepts, which makes their understanding more difficult. Although the latter explanatory model is more potent than the first one, because it is capable of explaining more cases of natural phenomena, e.g., motion involved in FS phenomena, concurrently, it is more complex, and hence, more difficult to understand, especially for younger students [4]. This is one of the reasons why primary and junior high schools adopt the density-based approach, in other words, the elimination of variables approach. In this case, teachers are faced with the difficult choice of introducing density using one of the following three ways: (a) the mathematical ratio of mass per volume, (b) the particle theory of matter, or (c) a visual representation that emphasizes the qualitative aspect of density [15]. However, the first case, i.e., a mathematical introduction of density, has been proven to be ineffective, due to the fact that students find it hard to differentiate between the concepts of mass, volume and density, or understand that density is an intensive quantity [15,34]. In addition, the second case, i.e., using the particle theory of matter, would most probably create misconceptions or reinforce students' prior ideas, such as the non-differentiation of density and denseness [43]. It has been shown that introducing density through this approach is rather abstract, which makes it very difficult for 9- to 12-year-old students to comprehend [44]. It would, thus, appear that the last case is the most appropriate because, in contrast to the first two, the use of a visual representation also provides opportunities for differentiating the focal concepts [26].

In the last few decades, teaching approaches have been developed in the framework of inquiry, emphasizing both the content of science and scientific practices [29,45,46]. As can be seen in Figure 1, there is a need for the emphasis to be on the understanding of the metaconcepts of scientific practices, e.g., the process of developing evidence-based conclusions through an experimental procedure, which, from an early stage, had been the focus of the elimination of variables approach. By this comment, we do not claim that inquiry-based teaching is not suitable or feasible to be implemented in the scientific approach of the analysis of equilibrium for FS phenomena. Rather, in order for the results of the eliminating of variables, experimental procedures to be understood and adopted by students, they need to be aware of the reasoning behind the scientific practices involved [15,46]. We claim that the TLS described in the next section is an example of effective density-based implementation, concurrently aiming at declarative knowledge (FS and density), procedural knowledge (control of variables strategy), and epistemological knowledge (nature of models) [18].

#### **4. Description of the TLS: Emphasis on the Content of FS and Its Didactic Transformation**

*4.1. The Design Principles of the TLS*

A case example of the elimination of variables approach for teaching and learning FS phenomena is presented with the five-unit TLS entitled "Density of materials in floating and sinking phenomena: Experimental procedures and modelling". The design principles of the TLS have been discussed in detail in former publications, e.g., [19]. For the sake of clarity, a brief reference to the six design principles of the TLS is made here, focusing on those regarding FS didactic transformation and relevant activities. More specifically, the design principles of the TLS were as follows [19]:


#### *4.2. The Didactic Transformation of the TLS*

FS phenomena are not included in the Greek primary school curriculum. However, the concept of density is introduced in fifth grade (10–11-year-olds) as a property of materials, with a limited number of examples, including one task that negotiates the sinking of a real ship. Although the curriculum proposes a guided discovery approach for negotiating phenomena and concepts, and one of the aims explicitly referred to is for students to understand the scientific method, the majority of teachers implement traditional deductive teaching-learning practices, followed by experimental demonstrations, whereas group work is sporadic [19].

As mentioned, the TLS was developed within the "Materials Science" project and is proposed as part of a broader curriculum for primary and lower-secondary level students (10- to 15-year-olds). The general objective of this broader curriculum is to restructure students' conceptual framework as regards the concepts of matter and material kind [33], including fluids. The elimination of variables rather than an analysis of equilibrium explanatory model for FS phenomena was adopted for reasons to do both with students' difficulties and the project's characteristics. Students' difficulty in understanding and effectively implementing an equilibrium mechanism, such as the buoyancy-based model to explain FS phenomena, has already been documented [4,15] and is analytically discussed in Section 3.

Furthermore, the emphasis of the project, which, on the one hand, is on the properties of materials properties, and the other, on inquiry-based teaching and learning, appears to be more compatible with the elimination of variables approach, that is, the density-based approach for teaching FS phenomena [15] (see Section 3). By highlighting the variables of these phenomena in order to derive a prediction rule that will determine "which" body will float, instead of negotiating the forces that are acting on an object when it is immersed in a liquid, we believe is easier for students in this age range to understand, and thus more conducive to the teaching/learning process.

It was decided to introduce the concept of density through the visual "dots-in-a-box" representation (Figure 2) as a property of materials [26]. As students would have already investigated the variables that affect the FS phenomenon, this visual representation shows several variables in only the one diagram. Obviously, the "dots-in-a-box" representation depicts the weight of the object by the number of dots, and the cube represents its volume or size, while its kind of material is now assigned to the conjunction of the weight and volume and not to any realistic representation of its external appearance. This enables students to easily make comparisons of the densities of the different materials in order to predict and explain FS phenomena, and if possible, to grasp the "heavy for its size" intermediate/precursor concept of density. In this way, the usual introduction of density using the mathematical ratio mass to volume, which has been shown to be difficult for students of these ages to grasp such relationships, has been bypassed [34].

**Figure 2.** The visual "dots-in-a-box" representation of the density of several materials, reprinted by permission from Springer Nature Customer Service Centre GmbH: Springer, Iterative design of Teaching-Learning Sequences, by D. Psillos and P. Kariotoglou, 2016 and by permission from [20], http://earthlab.uoi.gr/tel/index.php/themeselearn (accessed on 20 March 2021).

With a view to preventing or eliminating any misconceptions, such as "all hollow objects float" or "objects with air always float", as reported in the study by Yin et al. [29], we also thought it best to first introduce students to homogenous objects in the variable of solids. In the TLS, density is introduced within the context of floating or sinking of various everyday objects. First, students are introduced to homogeneous and then to composite objects, such as an iron cube and a ship, respectively, using causal relational reasoning, in other words, by comparing the density of the object and liquid, one is able to interpret and/or predict the FS of each object [31].

In sum, the didactic transformation of the TLS consists of two core choices (a) selecting the density-based approach to teach FS, and (b) qualitatively introducing the concept of density through the "dots-in-a-box" representation. Reference is made here to the other two subsidiary but significant aspects of the didactic transformation: (a) integrating scientific and technological knowledge into a context-based approach (fourth design principle), and (b) incorporating guided discovery experimentation into the explicit teaching of procedural and epistemological knowledge within an inquiry-based approach (fifth design principle). The teaching/learning environment for the experimentation activities was initially structured, which was gradually decreased with the aim of enabling students to become more autonomously involved [19,53].

More specifically, FS phenomena are introduced in the fourth design principle through a technological-problem scenario based on salvaging the Sea Diamond shipwreck [19]. The scenario, which runs throughout the entire TLS, has a dual role, on the one hand, it forms the familiarization phase, and on the other, it involves the following scheme: "technological problem", "scientific investigation", and "return to the problem" with the aim of increasing students' curiosity and motivation leading to the solution of the problem [49,50]. In contrast to the traditional approach, which proposes only the scientific investigation of the phenomena and the related concepts, our context-based approach through the technological-problem scenario increases students' interest and succeeds in involving them in the entire teaching/learning process [48].

The content of the TLS within the fifth design principle, which is the IBSE approach, includes elements of the inquiry method, i.e., aspects of the control of variables method as well as the nature and role of the models [18,19]. In other words, students are explicitly taught that to test if a variable influences a phenomenon, e.g., FS, then only this variable should differ, and all the other independent variables should be controlled [55]. In addition, students are explicitly taught aspects of nature and the role of scientific models, for instance, that models are not an exact representation of reality and that they are used to describe, predict, or explain a phenomenon [19,56]. Our hypothesis that procedural and epistemological knowledge would positively affect the understanding/interpretation of FS phenomena has been confirmed [18]. Thus, our claim that inquiry as a teaching goal constitutes part of the TLS's didactic transformation has been reinforced.

Summing up, the didactic transformation of the TLS described in this paper is based on (1) the elimination of variables approach to teaching FS phenomena, and (2) the qualitative introduction of the concept of density using the "dots-in-a-box" visual representation [26]. The explanatory model proposed for students to use in order to predict and explain FS phenomena for both solid and composite objects is based on the "dots-in-a-box" representation, in conjunction with the causal relational FS rule, which is, if an object's density is smaller than the liquid's density, then the object will float, and if an object's density is greater than the liquid's density, then the object will sink. We maintain that it is easier for young students (primary and junior high school) to grasp a qualitative representation of density, as a property of materials, rather than the scientific knowledge of the specific content, which is traditionally presented as a mathematical ratio of mass per volume, or even weight per volume [26,27]. Understanding the concept of density as "heavy for its size", thus perceiving density to be related simultaneously to both weight and volume, was an implicit teaching goal of the TLS. In this sense, we claim that by using the "dotsin-a-box" visual representation in their explanations for FS phenomena, as a property of materials and not of objects, students can differentiate the concepts of density and weight and consequently come closer to an intensive perception of the concept of density.

Every didactic transformation of content in science education is characterized by limitations related to (a) the kind of explanation of the focal phenomena and/or (b) the range of the phenomena explained when the model that is related to the transformed content is being used [12]. Therefore, the explanatory model of FS presented here (visual representation of density in conjunction with causal relational FS rule) has some limitations in comparison to other more abstract explanatory models (e.g., the analysis of equilibrium model for the explanation of FS and/or mathematical ratio for the representation of density). These are:


#### *4.3. The FS Content of the TLS*

The teaching and expected learning trajectory of the TLS implementation, along with a brief description and the activities in the five units, are presented here. The sequence of the activities is a fundamental element of the teaching design and implementation. In order to focus on the content that is directly relevant to FS phenomena, we do not include here any content to do with the control of variables strategy and models, which has, however, been described in previous publications [18,19]. The description follows the scheme of "main aim, content, and activities that students participated in" (Table 1). Students worked: (a) in groups of three to complete structured worksheets on both the real and simulated experiments that followed the POE (Predict–Observe–Explain) teaching strategy [57], and/or (b) in a whole class arrangement, following formative assessment activities [4].

In the first unit, our main aim was to provide students with a familiarization phase of FS through the technological problem of salvaging the shipwreck. Students participated in: (a) an introductory discussion about the variables that might influence floating and sinking, resulting in five independent variables: weight, material, and shape of the object, width of the container, type of liquid; and (b) a thorough discussion about the concept of a solid and homogeneous object, in contrast to a hollow object, in order to focus on solid objects. The teacher demonstrated an experiment in the POE approach to check if the first of the five independent variables, i.e., the weight of an object, influences the object's FS situation.


**Table 1.** Main aim, content and activities concerning FS in the five units of the TLS.

In the second unit, our main aim was for the students to understand that the FS of a solid and homogeneous object is influenced by the material of the object and the type of liquid. The POE activities that students participated in were characterized by the gradual decrease of scaffolding or the gradual increase of openness. In these activities, students were prompted to test the other four variables that might influence the FS of solid and homogeneous objects; we note here that the teacher has demonstrated the variable weight in the previous unit.

In the third unit, the main aim was the introduction of the "dots-in-a-box" representation of density and the use of this representation to predict and explain FS of solid and homogeneous objects, in conjunction with the causal relational FS rule, i.e., if an object's "dots-in-a-box" representation is smaller than the water's, then the object will float in water, and if an object's "dots-in-a-box" representation is greater, then it will sink in the water. Students: (a) searched and gathered information about the properties of several natural and artificial materials, such as glycerin and polyurethane, and the ways they can be used, (b) negotiated with cubes of the same volume but different material and were assigned with a task prompting them to express the "heavier-lighter" material relationship (Figure 3), and (c) completed a task in a simulated environment, using a balance in order to put cubes of the same volume but different materials in the order of heavier to lighter. The sequence of densities of the materials in Figure 2 resulted from this sequence of tasks.

**Figure 3.** One of the students' proposals to represent the "heavier-lighter" relationship of cubes of different materials, noted on the relevant worksheet.

In the 4th unit, the main aim was to generalize the causal relational FS rule to explain FS in any liquid and for both solid and composite objects. First, students participated

in real and simulated experiments of several homogeneous objects in glycerin instead of water, emphasizing once again the role of the type of liquid in the explanatory model of FS phenomena (Figure 4a). In addition, understanding that the density of a two-material composite object lies between the densities of these two materials was crucial, so students participated in real experiments of composite objects in water, emphasizing the role of the average density of an object in predicting and/or interpreting its FS. In this unit, the phrase "dots-in-a-box" was replaced with "density" to refer to the property of materials.

**Figure 4.** (**a**) A simulated experiment using glycerin instead of water, reprinted by permission from Springer Nature Customer Service Centre GmbH: Springer, Iterative design of Teaching-Learning Sequences, by D. Psillos and P. Kariotoglou, 2016; (**b**) simulated environment to investigate FS of the Sea Diamond shipwreck.

> In the 5th unit, the main aim was to find the best solution for lifting the shipwreck, if possible, by implementing the generalized causal relational FS rule. The students had the opportunity to work in groups both in real and in a simulated environment (Figure 4b) so as to investigate the FS of the Sea Diamond shipwreck, studying the effect of excess water in the hold of the ship and discuss how it could be salvaged.

#### **5. Selected Results**

Here, a short description of the research method, and selected results on the learning of FS, of the implementation of the revised TLS are presented to show its effectiveness and its didactic transformation.

#### *5.1. Method*

The TLS described in this essay was developed, implemented and evaluated twice in an iterative evolution manner [19]. The initial and the revised versions were implemented on twelve and forty-one 5th graders, respectively, in a real-classroom educational context [20–22].

Data were collected using several sources (questionnaires, interviews, researchers' notes, etc.). For the sake of brevity, we focus on six Tasks of the questionnaire (see Appendix A) in order to evaluate the impact of the TLS implementation on students' FS explanations. The students answered the same questionnaire before, immediately after the TLS intervention, and seven months later (Pre, Post, and Delayed Post, respectively). More specifically, Tasks 1 and 2 were answered all three times, whereas Tasks 3–6 were answered in the Post and Delayed Post Tests. These tasks were not included in the Pre-questionnaire, either because students did not know the 'dots-in-a-box' representation before the intervention (Tasks 3 and 4) or because the questions were too complex to be attempted before the intervention (Tasks 5 and 6). While the results for Tasks 1–4 have been published in previous papers [20,21], the findings for Tasks 5 and 6 are presented here for the first time.

#### *5.2. Results*

The analysis results of students' responses for the six Tasks (Means and Standard Deviations) are presented in Table 2.


**Table 2.** Means and Standard Deviations of the students' responses in the six Tasks (*n* = 41).

Tasks 1 and 2 examined students' explanations about FS phenomena, whereas the 'dots-in-a-box' visual representation was not given. In both Tasks, the weight or the size of the object is the characteristic that could mislead students' responses to the alternative idea that "an object sinks because it is heavy" or "because it is big". The responses were classified thus: 0 for no or irrelevant answers, 1 for reference to the object's weight, 2 for reference to the material the object consists of, and 3 for reference to the causal relational FS rule, i.e., comparing the densities of the object and the liquid to predict and/or explain the object's FS. The study findings showed that there was a statistically significant improvement after the TLS intervention for both Tasks 1 and 2 (*z* = 3.446, *p* < 0.001 and *z* = 3.801, *p* < 0.001, respectively), which was retained seven months later. More specifically, for Task 1, the mean was 1.10 for the Pre-Test, which went up to a high 1.78 in the Post-Test, and which was maintained with a slight decrease seven months later in the Delayed Post-Test (1.68). The results were similar for Task 2, where the mean in the Pre-Test was only 0.78, which rose to 1.39 immediately after the TLS intervention in the Post-Test and was maintained almost at the same level (1.37) in the Delayed Post-Test.

Tasks 3 and 4 examined students' explanations about FS in a simulated environment, while the "dots-in-a-box" representation was also given. Students' responses were classified as 0 for causal linear reasoning, i.e., focusing only on one characteristic of the objects, and 1 for causal relational reasoning, i.e., a comparison of the densities of the object and the liquid. The results showed that most of the students were able to use the "dots-in-abox" representation to successfully apply causal relational reasoning in their responses, immediately after the intervention and seven months later (Table 2). More specifically, for Tasks 3 and 4, the means for the Post-Test were 0.63 and 0.73, respectively, which decreased slightly in the Delayed Post-Test to 0.56 and 0.54, respectively. In addition, it appears that most of the students who were able to apply causal relational reasoning in Task 4 could also understand differences in float levels in relation to material density [20], a topic that had not been covered in the intervention.

The last two Tasks 5 and 6 examined whether students could effectively determine the position of an object, in relation to the surface of a liquid, in order to establish its density. In Task 5, in order for students to correctly decide the relationship between the densities of the two objects by applying the causal relational FS rule, they would have to disregard the size of the two objects, which could have been misleading. Students that did not, even qualitatively, differentiate between weight and density would intuitively think that the bigger object had a greater density. Students' responses were classified as 1 for using the causal relational FS rule and 0 in all other cases. The results in Table 2 show that a large number of students were able to successfully answer Task 5, with a mean of 0.44 in both the Post- and delayed Post-Tests, thus indicating that they, at least, qualitatively differentiated between weight and density. For Task 6, in order to successfully apply the causal relational

FS rule and come to the conclusion that the density of the object is equal to the density of the liquid, students had first to recognize that the object was in a state between floating and sinking, and therefore suspended, which meant that it remained in the liquid at the same location at rest. Responses were classified as 1 for recognizing that the object was suspended and for successfully applying the causal relational FS rule, and 0 in all other cases. The results in Table 2 show that a large number of students were able to successfully answer this question, with means of 0.51 and 0.46 for the Post- and Delayed Post-Tests, respectively. This finding indicates that the students were able to apply the causal relational FS rule in a new situation, that of the suspension of an object in a liquid, which had not been covered in the TLS intervention.

#### **6. Discussion and Conclusions**

In this paper, we describe and examine certain elements in the teaching of floating and sinking, with focus on the didactic transformation of content, which was part of a long-term developmental study. The didactic transformation process is an essential aspect in each teaching effort and especially in the design and development of a TLS.

In our research on the teaching of floating/sinking (FS) phenomena to 5th grade Primary school students, we adopted the elimination of variables approach (density-based) rather than an analysis of equilibrium explanatory model (buoyancy-based). The reasons for this decision were: (1) students find it difficult to understand and effectively implement the buoyancy-based model to explain FS, in contrast to density which is an easier concept to comprehend; (2) The emphasis of our developmental study project was on materials and their properties; and (3) Inquiry-based teaching and learning is more compatible with the elimination of variables approach, for the teaching of FS phenomena, as well as being more in line with the current Greek school curriculum.

The concept of density was introduced to students qualitatively through a visual representation called "dots-in-a-box". The qualitative method was chosen over (a) the mathematical ratio of mass per volume and (b) the particle theory of matter. The reasons for this choice are as follows: (1) the visual representation makes it easier for students to grasp this scientific concept, whereas the mathematical introduction of density has several times been shown to be ineffective, at least with primary and junior high school students [34]; (2) misconceptions would most likely arise with the particle theory of matter; and (3) introducing density with a visual representation also provides students with opportunities for differentiating the concepts involved in the interpretations of FS phenomena, e.g., density, weight and volume.

The study findings strongly suggest that our TLS and the didactic transformation of content that was developed had a significant level of success in the teaching/learning of FS to young students, which seems to have been maintained seven months later. Most of the students adopted explanations that were compatible with scientific ones and were able to overcome their prior alternative ideas, such as "heavy objects sink and light objects float". In addition, when given the "dots-in-a-box" visual representation of density, the students successfully implemented the causal relational FS rule, i.e., comparing the densities of the objects and liquid to predict and/or explain the FS state of objects. Finally, several students were able to implement the density-based model, which they had been taught, to explain situations that had not been covered in the TLS intervention, such as successfully predicting the floating level of objects made of different materials and applying the causal relational FS rule for objects suspended in a liquid. Both cases are considered difficult for students in this age range (10–15 y-o) to comprehend and explain [1].

We consider that this work is an original sample of a developmental research description, in the sense of the Design-Based Research approach in science education [8,23,24], and consequently an example of effective good teaching practice, that can help teachers to elaborate on their teaching and inspire innovative treatment of the topic of floating and sinking in science curricula for primary school physics. The teaching and learning intervention for floating/sinking phenomena, which is in itself a difficult conceptual science topic, to this

young target population, was successful, we strongly believe, due to the contribution of the didactic transformation of content.

**Author Contributions:** The authors contributed equally to the preparation and writing of the present paper. All authors have read and agreed to the published version of the manuscript.

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

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **Appendix A**

The questionnaire tasks

*Task 1*

On a big ship, among other objects, you can find an anchor. Does it float or sink if we drop it into the sea? Justify your answer.

... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

The anchor: floats sinks I do not know

Because:

Costas drops a small piece of a particular material into a container filled with water, and he observes that it floats. Afterwards, Irene drops a bigger piece of the same material into the same container. In your opinion, at which point will the big piece stop moving? Circle which number: 1, 2 or 3 in the diagram you think represents the final position of the two bodies that Costas and Irene dropped into the container. Justify your choice. *Task 3*

You are given two objects A and B, and a container which contains a liquid. The densities of the two objects and that of the liquid are given with the "dots-in-a-box" representation, as you can see in the gray box. If you drop objects A and B into the container with the liquid, what will their final position be? Draw objects A and B in their final position in the liquid. Justify your answer:

... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

The densities of the two objects and that of the liquid are given with the "dots-in-a-box" representation, as you can see in the gray box. We drop objects A and B into the liquid. Circle which number: 1, 2 or 3 in the diagrams best represents the final positions of the two objects after we have dropped them into the liquid.

Justify your choice:

... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

*Task 5*

We drop the two Objects (A) and (B) in the liquid (C). Object (A) floats in the liquid, whereas Object (B) sinks in the liquid. Decide if the following sentence is correct or incorrect:

Object (A) has a greater density than Object (B).

It is right It is wrong I don't know

Justify your choice:

... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

*Task 6*

While Georgia, Petroula, and Sophia were playing with some toys, one of the toys accidentally fell in the container with the liquid that you can see in the picture. The girls noticed that the toy did not float up towards the surface of the liquid, nor did it sink to the bottom of the container. They wondered what the density of the toy could be, but they disagreed in their opinions:

Georgia says that this object has a greater density than the liquid.

Petroula believes that the object has a lower density than the liquid.

Sophia says that the object has the same density as the liquid.

Which of the girls do you agree with? With:

Georgia Petroula Sophia I don't know

Justify your choice:

#### ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

#### **References**


### *Article* **Open Water Flume for Fluid Mechanics Lab**

**Rachmadian Wulandana**

Mechanical Engineering Program, SUNY New Paltz, New Paltz, NY 12561, USA; wulandar@newpaltz.edu

**Abstract:** Open water flume tanks with closed-loop circulation driven by centrifugal pumps are essential for hydro experimentation in academic settings as well as research centers. The device is also attractive due to its versatility and easy-to-maintain characteristics. Nevertheless, commercial open flume systems can be expensive and become less prioritized in engineering schools. This paper describes the design and fabrication of an affordable, medium-size water flume tank, suitable for education purposes. The central piece of the system is a transparent observation chamber where fluid experiments are typically conducted and observed. The expected maximum average water speed in the observation chamber of about 60 cm per second was achieved by the inclusion of a 3 hp centrifugal pump. The size and capacity of the current design were constrained by space limitation and available funds. The educational facility was assigned as a two-semester multi-disciplinary capstone senior design project incorporating students and faculty of mechanical, electrical, and computer engineering programs in our campus. The design process provides a training platform for skills in the area of Computer Aided Designs (CAD), Finite Element Analysis (FEA), Computational Fluid Dynamics (CFD), manufacturing, and experimentation. The multi-disciplinary project has contributed to the improvement of soft skills, such as time management, team working, and professional presentation, of the team members. The total material cost of the facility was less than USD 6000, which includes the pump and its variable frequency driver. The project was made possible due to the generous sponsor of the Vibration Institute.

**Keywords:** open water tank; education; fluid mechanics

#### **1. Introduction**

Open water flumes provide effective hands-on multidisciplinary learning tools for students matriculated in engineering programs, as well as in the physical education program [1]. Senior students participating in the design and fabrication process of the flume obtained the opportunity to reinforce fundamental engineering concepts and to master valuable technical skills [2–4]. The reported open flume was motivated by the need of a test chamber for an ongoing research on bladeless turbines, as well as the desire to have such device in the fluid mechanics laboratory. The current design was inspired by a small scale commercial water flow tank used in the investigation of vortex-induced autorotation and oscillation of straight cylinders [5,6], as well as the energy potential from such vibration modes of symmetric geometries [7,8]. Open channel tanks are paramount for various hydrodynamic research areas, such as designs of hydrokinetic energy harvesters, investigation of flow characteristics in the presence of obstacles, drag and lift of objects exposed to fluid flows, etc. The vortex-induced vibration (VIV) of objects exposed to flow has emerged as attractive potential sources of renewable energy. Additionally, the versatility of open channel flumes has allowed plethora of vibration modes to be experimented. Sun et al. used a 40 cm wide and long rectangular channel to convert vortex-induced vibration and galloping of blunt bodies into electricity via piezoelectric strips [9]. In a small scale water tank, Cao et al. attempted to amplify the vibration amplitudes of piezoelectric strips by means of magnet fields [10,11]. Arionfard and Nishi utilized a small (30 cm × 1 m) water channel to investigate power harnessing potential of pivoted cylinder exposed to highly turbulent flow [12]. The utilization of much larger size channels allows modeling of

**Citation:** Wulandana, R. Open Water Flume for Fluid Mechanics Lab. *Fluids* **2021**, *6*, 242. https://doi.org/ 10.3390/fluids6070242

Academic Editor: Mehrdad Massoudi

Received: 31 May 2021 Accepted: 26 June 2021 Published: 3 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

near-realistic flow and full-scale energy harvesters. Rostami and Fernandes studied the energy harvesting from the torsional galloping, fluttering, and autorotation modes of flat and S-shape plates using a 1.4 m, wide 22 m long water channel [13–16]. Furthermore, a 1 m wide water flume tank is used in the development of VIVACE (Vortex-Induced Vibration of Aquatic Clean Energy), which exploits the transverse-to-the-flow (lateral) oscillation mode of VIV at large range of Reynolds numbers using magnets attached to cylindrical bars as harvesters [17–19]. An open channel with similar width but five times longer (about 42 m long) was recently used in the attempts to utilize the galloping modes of triangular cylindrical bars at much higher velocity [20,21]. Generally, interests in the designs of innovative hydrokinetic energy turbines have been well facilitated by open water flumes. The facility is adaptable for various types of turbines and its blades, as well as a range of measurement needs. Barber et al. studied the power and thrust improvement of axial turbines after the utilization of carefully designed adaptive pitch composite blades [22]. A small commercial flow loop water tank was used in the investigation of wake behind the propellers of wind turbines [23]. Open water flume experimentations of vertical axis hydrokinetic generators, such as Savonius and Darrieus turbines, have been intensified lately. These turbines are attractive renewable energy harvesters due to their simple designs, low-cost manufacturing, easy maintenance, and independence of flow direction. Talukdar et al. demonstrated that the performance of two-bladed semi-circular Savonius turbines is better than two-bladed elliptical and three-bladed semi-circular designs [24]. On the other hand, Sarma et al. pointed out that the three-bladed Savonius model performs better in water than in air flow [25]. It is interesting to note that placing an upstream obstacle can improve the performance of a modified two-bladed Savonius turbines [26] and Darrieus turbines [27]. Attempts to assemble and test the two turbines into one single hybrid unit have shown promising outcomes [28,29]. In the hydrology area, large size tanks facilitate investigations of particle sedimentation carried by water flow [30], improvement of hyporheic zones due to riverbed restoration [31], and effects of flood waves on river banks [32] among many other studies. The open channel allows examination of flow qualities, such as head loss and turbulence, by placing out obstacles in its test section. Investigation of the head loss due to the presence of submerged baffle-posts [33,34] and examination of water turbulence due to rib roughness [35] are examples of important studies to be carried out for better irrigation systems. Examinations on the interaction between vegetation and water flows shed light on the life quality of river. Here, the vegetation under investigation is carefully arranged along the base of the test chamber and its effects on water waves are studied using camera [36]. Lastly, we would like to mention flume experimentations to test the drag and lift of 3D-printed models of swimmers [1] that show the capability of water flume to accommodate large range of multidisciplinary applications. The transparent chamber of the flume certainly allows flow measurement by either direct observation [33,34] or flow visualization technology such as dye injection [31], hydrogen bubble [37], particle image velocimetry (PIV) [23,37,38], planar laser-induced fluorescence (PLIF) [39], stereoscopic digital particle image velocimetry (SDPIV) [40], laser doppler velocimetry (LDV), and acoustic doppler velocimetry (ADV) [30] methods. It is obvious that the versatility and easy-to-maintain characteristic of the open flume systems makes the device very attractive to be included in fluid mechanics and hydrology labs.

The presented work discusses the design and fabrication of a small size (15 cm wide and 75 cm long test chamber) open water flume with closed flow loop driven by a 3 hp centrifugal pump. The designs must consider limited funds, available space in the lab, and mobility requirement. The overall size of the system, the materials for the structural frame and the chambers, as well as the pump selection, can be determined by the flow requirement in the observation chamber and experiments that would be conducted. In this project, the maximum fund was set to be USD 10,000 and the overall size is governed by the size of the freight elevator in the school. The overall length of the system must be kept to about 2.5 m or shorter, so that it can fit in the school's freight elevator. Reported in this paper are the pump selection process, structural analysis of the supporting frame, computational fluid dynamic analysis of the flume, manufacturing process, and conducted experimentations on the flow visualization using dye injection and investigation of drag forces on submerged objects. A brief discussion on the budget is presented before the paper is closed with a discussion and conclusion sections.

#### **2. Basic Design and Pump Selection**

This water flow tank is designed to provide straight uniform flow with an average speed of as much as 60 cm/s through its transparent observation chamber. The crosssectional dimension of the observation chamber is planned to be at least 15 <sup>×</sup> 15 cm2. The length of the observation chamber was designed to be 60 cm. This length is needed to provide enough observation room for vortex trail and wakes behind objects exposed to water flow. The final length of the test chamber is 75 cm. The selection of the maximum speed and chamber dimension was based on a commercial water flume used in experiments on vortex-induced oscillation of cylinders [41]. An upstream manifold in a shape of converging chamber of about 44 cm is added to provide room for the water to reduce its turbulence and complex characteristics prior to enter the observation chamber. The height of the observation chamber was designed to be ~1.2 m, a little less than the average height of human eyes of ~1.4 m. The top side of the observation chamber is expected to be open to allow direct physical access to the flow and, more importantly, easy placement of objects exposed to the flow. Tight covering of the top side to achieve high speed of water flow is possible, but the pressure increase must carefully be calculated. The flow tank was designed to be mobile, so that it can be easily moved when needs arise. The overall length and width of the flow tank therefore is constrained to size of the school's freight elevator of approximately 2.16 by 2.44 m2. This size put limitation on the total length of the equipment to be about 2.5 m. The overall budget of this device is set to be less than USD 10,000, based on the maximum fund by the sponsor for this project: Vibration Institute. As the project is scheduled for a senior design project, the design and construction must be finished in two semesters.

The design process begins with a preliminary calculation for the pump specification, which would determine the pipe diameters, the dimension of supporting frames, and required power outlet. Detail calculation would require flowrates, height of the observation chamber, lengths, diameters, and materials of various pipes involved in the designs, as well as connections and pipe bends. In the absence of many of these parameters, the needed pump power can be calculated using the required flowrates and maximum elevation that needs to be overcome by the water. Other parameters need to be estimated. Figure 1 panel (a) shows an example of simple conceptual sketch of the water flow system drawn by students involved in this project. The flowrate can be determined from the demanded speed *V* of ~0.6 m/s and cross section area of the observation chamber of *<sup>A</sup>* <sup>=</sup> 0.15 <sup>×</sup> 0.15 <sup>=</sup> 0.0225 <sup>m</sup>2, through the equation for flowrate *<sup>Q</sup>* <sup>=</sup> *<sup>V</sup>* <sup>×</sup> *<sup>A</sup>* <sup>=</sup> 0.0135 m3/s, or about 12,800 gph or 214 gpm. Figure 1 panel (b) depicts a simple open flow diagram to represent the design of the water tank. The pump must overcome the height of the observation chamber and resistance by all pipes and chambers (major losses) and various connections (minor losses) along the loop. Note also that the effects of gravity on the returning channel from the observation chamber is not included in the calculation, as the power discount may not be significant. The flow diagram shown in Figure 1 guides us in determining the Bernoulli equation to be used.

**Figure 1.** (**a**) Sketch of 1st design of the flow tank loop given by students and (**b**) the basic flow diagram for the calculation of pump head needed to determine the pump power.

The total head *hp* for the calculation of pump power can be obtained from the Bernoulli equation easily found in fluid mechanics text books [42,43];

$$h\_p = \Delta h + \Sigma \, f \frac{L}{D} \frac{\nabla^2}{2g} + \Sigma \, K \frac{\nabla^2}{2g} \tag{1}$$

The Δ*h* represents the elevation difference between two end points in the flow diagram. In our case, we would take this as the largest height difference available in the water tank system, which is the height of the observation chamber. The second and third terms in the equation are known as the major and minor losses, losses that are caused by the wall friction along the pipes and channels in the system, and losses that are caused by any other obstacles in the system, respectively. For the major loss, the Darcy's friction factor *f* is known to depend on the Reynolds numbers of the flow and roughness of the pipe. The *L* and *D* are length and diameter of pipe, respectively. In the third term, the constant *K* represents resistance factors for various obstacles along the loop, such as bends, junctions, entrances, etc. The Reynolds number is defined as

$$Re = \frac{\rho \nabla D}{\mu} \tag{2}$$

where *ρ* and *μ* are the water density and absolute viscosity, respectively. In calculating the major and minor losses, the average speed *V* of the water should be taken at each section along the pipe where the loss is calculated, and the *D* represents hydraulic diameter of that particular section. As the detail geometry of the flume is not yet known, we assume the contribution of these losses to be simply double the elevation difference, hence the total head is *hp* = 1.2 + 2 × 1.2 = 3.6 m. It will be shown later that this estimation is larger than an alternative estimation based on possible major and minor losses in the flume. The required pump power can be calculated using [42,43]

$$\mathcal{W}\_p = \frac{\rho \text{ g } \mathcal{Q} \, h\_p}{\eta} \tag{3}$$

In estimating the needed pump power, we took the density (*ρ*) and absolute viscosity (*μ*) of water at atmospheric pressure and 20 ◦C to be 998 kg/m3 and 0.001 Pa.s [43,44], respectively, while the gravitational constant *g* = 9.81 m/s2. Not knowing the pump to be purchased yet, we estimated the pump efficiency *η* from commercial pump examples presented in textbooks [43,44]. Based on these examples, for the required flow rate of 214 gpm, the efficiency was found to be about *η* ≈ 60%. Using these data, we can estimate the required pump power to be only *Wp* = 793 W or about 1.06 hp. Note that the power

calculation here is based on high flow rate but very low head loss, due to the very short pipe and minimal height involved in the current design. Consequently, such a pump, with such combination of low pump power and high flow rate, may not be available in the market. Here, students need to realize that the market availability of the pump puts a constraint in the design.

Alternatively, the major and minor losses could be estimated based on the possible materials and components to be used in the designs. This detail calculation process provides valuable realistic exercise for the students involved in this project. Table 1 shows list of pipes and channels, as well as estimated major head losses. The mean velocity for each section is estimated from the flowrate of 0.0135 m3/s based on the required flow speed of 60 cm/s. The major loss is defined as *h* = *f <sup>L</sup> D V*2 <sup>2</sup>*<sup>g</sup>* and the friction factor for each section is determined from the Moody chart based on the associated Reynolds number and given roughness factors of the pipes [42–45]. Due to the short pipes involved in this design, this calculation renders net head much less than 1 m, well below the estimated head mentioned before. As expected, the required head is dominated by the height of the water column that must be overcome.

**Table 1.** List of friction and major head loss of pipes and channels involved in the designs.


Similarly, the minor losses can be estimated by considering all possible pipe components that are to be included in the system and evaluate its contribution. Table 2 shows the list of possible components to be installed in the system and its corresponding "*K*" values obtained from a textbook [42,45]. The total *K* value is about 4, and this results in minor head loss of about 0.2 m, after assuming velocity of 1 m/s. Combined with the major head loss calculated above, the total would still be less significant compared to the head by elevation difference.

**Table 2.** In this table, we list of possible minor loss components to be included in the design.


Three (3) pumps for minimum 190 gpm flow rate are selected from McMaster– Carr [46]; (A) High-Efficiency Circulation Pumps for Water, Coolants, and Oil, (B) Harsh-Environment Self-Priming Circulation Pumps for Water and Coolants, and (C) High-Flow Inline Circulation Pumps for Water. Shown in Table 3 are features of the three pumps considered for the system. The price range is approximately USD 1000 to 2000 (2018 price). The specification of these pumps shows that the flow rates and power requirement slightly above our design requirement and the prices are within the budget range.


**Table 3.** List of potential centrifugal pumps available in McMaster–Carr and its working fluid, maximum flow rate, power requirement, maximum head, and unit price.

> In selecting the pump, students may be asked to setup a simple decision matrix. The matrix allows a group of students to quantify their opinions on factors that weigh the buying decision such as price, capacity, pump power, and type of fluids that can be handled by the pump. In using the matrix, each student in the design team gives score between 1 to 3 (as there are three pump candidates), indicating the level of preference, for each category. For example, a student who prefers to buy pump B as they think that the price is the most reasonable, not necessarily the cheapest, should score 3 for pump B under the category of "Price". If the student considers that they prefer another pump for the power, then they can put a score of either 2 or 1 under the "Power" category for pump B. Students involved in the pump selection process collect their matrix, and the pump with highest total score should be selected. Table 4 shows an example of such decision matrix performed by one student. This example shows that the student preferred to purchase pump B over the other two pumps.

> **Table 4.** An example of decision matrix for pump selection process. Each student involved in the project may participate in the selection process by filling out this table.


The selected 3"-self-priming pump is manufactured by AMT, a Gorman–Rupp Company based in Mansfield, OH, USA [47]. This is a 3-phase pump that draws 10 Amps at 208–230 V. The stainless-steel impeller is encased in cast iron and capable of delivering 375 gpm maximum capacity of water at its 3450 rpm. The pump is selected due to its availability, capability in delivering 200 gpm, and affordability. Testing of the pump after the complete assembly indicates that the flow speed reaches 60 cm/s at the maximum 50 Hz input.

#### **3. Static and Dynamic Analysis of Supporting Frame**

Structural analysis is paramount to determine the integrity of the flow tank and its supporting frame. Figure 2 panel (a) shows the final Computer Aided Design (CAD) drawing prepared by students and its major components. The final product of the water flow tank is shown on panel (b). The final total length of the water tank system is 2.51 m, which includes the 75 cm long observation chamber (1) and inlet manifold into the observation chamber (5) with enough entry length. In the discussion below, only the structural analysis of the supporting frame is presented.

**Figure 2.** Left panel (**a**) shows the drawing of the water tank and its major components; 1. observation chamber, 2. centrifugal pump, 3. variable frequency drive or pump controller, 4. return pipe to the pump, 5. inlet manifold into the observation chamber, 6. holding tank, 7. structural frame, 8. wheels, and 9. ball valve. On the right, panel (**b**) shows the final product of the water tank supported by frames and wheels.

The material chosen for the supporting frame is ASTM A513, which has a yield strength of 220.6 MPa (32,000 psi). Considering a safety factor of 1.3, the maximum yield strength for the system is *σ<sup>y</sup>* = 169.72 MPa (24,615 psi). The beam profile used is 5.08 <sup>×</sup> 5.08 cm<sup>2</sup> hollow square tubing with a thickness of 2.11 mm. The static analysis involves calculations of dead loads acting on the system, which includes the weight of the water and the tank. Table 5 lists all possible static load items and its estimated amount. The list was prepared by students as part of the full static and dynamic structural analysis of the supporting frame. Several possible designs were analyzed and results from one design is presented here.


**Table 5.** The table shows main loads to be considered in the design of the supporting frame.

The maximum deflection of 0.0254 mm was discovered on one of the cross beams directly supporting the observation chamber. This value is considered very minimal, and it indicates that the selected channel is sufficiently strong. In this model, the frame is assumed to be supported by pins on its four bottom corners. Static loads, representing items listed in the Table 5 are applied as point loads on several locations on top beams. The water weight is distributed throughout several points on the frame. Figure 3 displays the deformation of the frame in the full 3D static analysis. As is expected, the bottom beams show minimum deformation due to the supports, while the top beams show large deformation due to the application of load.

**Figure 3.** The deformed shape of the supporting frame due to the dead loads indicate maximum deflection of about 0.001 inch occurs on the cross beam that supports the observation chamber. The 3D static analysis in ANSYS is performed by members of capstone senior design team.

Additionally, a dynamic modal analysis was performed using ANSYS (ANSYS Inc., Canonsburg, PA, USA) to check the mode shapes and natural frequencies of the supporting system. Outcomes from the modal analysis indicate possible resonance and amplification of deformation when the pump's operating frequency is in the vicinity of the natural frequency of the support system. Due to the simple vibration load that is being considered, the modal analysis was set to find only the first six (6) natural frequencies. The resulting natural frequencies of the frame are presented in Table 6. These natural frequencies were used in harmonic response analysis to determine the maximum deflection of the frame. The maximum deformation of 5.51 mm occurs for Mode 5 corresponding to the frequency of 73 Hz. Other frequencies result in maximum deformation of less than 2.54 mm. The first three (3) modes correspond to the operating frequencies of the pump of 0 to 60 Hz. Nevertheless, due to the rigidity of the structure, the amplification was found to be very minimal. The outcome suggests that the structure and material selection for this flow tank is sufficiently strong to sustain the vibration load by the pump. The study also shows that only the first 3 (or at most 4) modes need to be considered as the maximum operating frequency of the pump is only 60 Hz. Table 6 shows the six (6) first natural frequencies of the supporting frame. The maximum deformation from the harmonic analysis, based on the given static load, is also shown in this table. A maximum deformation of 5.08 mm can potentially occur for Mode 5, corresponding to ~73 Hz. Other modes show maximum deformation of less than 2.54 mm.

**Table 6.** The first six natural frequencies of the supporting frames are listed in this table, along with the maximum deformation associated with the harmonic analysis.


The first, second, and third modes represent predominantly lateral deformation of the top beams due to the rigid support of the bottom beam. The lateral deformation can occur harmoniously between the parallel top beams (Modes 1 and 3), but it can also occur non-harmoniously as in Mode 2. Mode 3 shows lateral deformation in the long direction of the frame. While the dynamic analysis shows sufficient stiffness and unlikeliness for large amplification, rubber footings were placed underneath the pump's platform to damp the vibration.

#### **4. Computational Fluid Dynamics**

Computational fluid dynamics (CFD) analysis suggests useful information for the design of various channels of the flume. Nevertheless, the CFD simulation and analysis provide hands-on exercise materials that allow students to better understand the dynamics of flow phenomena through visualization and parametric studies [48]. Possible analysis to be performed ranges from two-dimensional to full three-dimensional study of either unsteady or steady water flow. Examples of topics to be studied include studying the flow characteristics in the observation chamber, pressure losses across the flume, undesired circulation zones, areas of turbulence, extreme stresses caused by fluid flow, etc. The computation domains of such studies can be either developed from scratch or imported from the structural mechanics CAD designs. In this section, results from a full 3D steady state analysis of water flows in a semi-loop model of the tank are briefly presented. The study is performed using the commercial package COMSOL Multiphysics 5.5 with CFD module (COMSOL Inc., Burlington, MA, USA). Most of the computations are performed using HP ProBook 640 G2 Notebook PC (Hewlett-Packard Company, Palo Alto, CA, USA). operated using 64-bit Windows 10 Enterprise version 1909 with a total capacity of ~465 GB. The computer has an installed RAM of 16 GB and is equipped with Intel® Core™ i7-6600U CPU @ 2.60 GHz. One model is run using a desktop HP Compaq Elite 8300 CMT computer equipped with Intel® Core TM i7-3770 CPU at 3.40 GHZ and 16.0 GB RAM.

The computation domains can be differentiated into several parts, as shown in Figure 4 panel (a). The pipe inlet domain represents the vertical cylindrical thick PVC pipe channel that delivers water from the pump's exit into the main part of the water tank. This channel is 8 cm in internal diameter and the length is 40 cm. The pipe is extended into the diffusing chamber before it splits into two horizontal branches with the same internal diameter. The resulting T-pipe allows the water to be distributed in the diffusing chamber. The T-pipe also prevents the water to "shoot up" and put high pressure on the ceiling of the diffusing chamber. The next important domain is the converging chamber that functions to collect the water and channel it into the observation chamber. A flow straightener could be included in the converging section, however this part is not modeled due to its large amount of mesh. This simplification is justified by the computation results that indicate straight uniform flow in the observation chamber (Figure 6). The observation chamber is a 60 cm long prismatic pipe with a square cross section of 15 <sup>×</sup> 15 cm2. At the end of the chamber, the water is collected in a rectangular drainage tank of 60 <sup>×</sup> <sup>30</sup> <sup>×</sup> <sup>40</sup> cm3. The CFD model does not include a horizontal transparent plastic pipe that returns the water back from the drainage chamber into the centrifugal pump.

**Figure 4.** Panel (**a**) shows the computational domains involved in the analysis and panel (**b**) shows the non-uniform mesh generation implemented for the analysis.

Several simplifications are implemented in the CFD model. For instance, the actual drainage chamber is converging downward, but in the model only the top 40 cm section is included. It is assumed that the inclusion of this part is sufficient to recreate necessary backflow effects that occur in the observation chamber. The converging chamber attached to the distal face of the observation chamber is made to be plain horizontal, parallel to the observation chamber. The actual configuration of this curved chamber includes a slight decline in the vertical direction into the entrance of the observation chamber. Lastly, in the CFD, the observation chamber is assumed to be a closed rectangular pipe. The actual observation chamber is open on its top side to allow direct access to the water.

The inlet boundary condition is located at the free end of the inlet pipe that represents the conduit between the centrifugal pump and water tank body. Not knowing the exact velocity profile produced by the centrifugal pump, a uniformly distributed velocity profile is assumed at the inlet. The exit boundary condition is located at the bottom of the drainage chamber and a zero-pressure outlet is assumed here. The actual pressure level required by the centrifugal pump to produce the flow can be obtained by adding known hydrostatic and dynamic pressures occurred at the exit level. The internal wall is assigned to the wall of the T-pipe domain, as this pipe is located inside the diffusion chamber. The physical T-pipe has a non-zero pipe thickness that has been ignored in this CFD modeling. This important feature certainly leads to a convenience and efficient mesh generation. The remaining walls are assumed to have zero velocity. Figure 4 shows the complete model, domains, boundary conditions, and partial mesh generation that results in non-uniform mesh size throughout the domains.

CFD analysis requires a delicate balance between demands for accuracy of the outcomes versus the computational cost that are the direct results of both mesh generation and complexity of the selected flow model. Table 7 lists CFD models studied in this work and their mesh generation scenarios. Models As and Bs are all meshed using "Global Mesh" method, where the same mesh size or scheme is applied uniformly on all domains of the model. However, models As are computed by means of "Laminar" flow model, while models Bs are executed using "Turbulent" flow model. The C, D, E, and F models are discretized using "Partial Refinement" method and are executed using "Turbulent" flow model. Effects of the two flow models will be explained later in this section. All models are computed using the HP ProBook laptop, except model D that is executed using the desktop computer. The Global Mesh method is straightforward, but for three-dimensional models this method may result in astronomical element numbers and computation time with possible less-than-sufficient accuracy at places of interest. In the current work, we later show that partial mesh refinement on places of interest should be practiced as a choice of practice in CFD modeling, particularly for three-dimensional modeling. In our work, the "Extremely Coarse" mesh option of COMSOL Multiphysics (COMSOL Inc., Burlington, MA, USA) for our model results in 21,011 elements. Increasing the accuracy by selecting the "Normal" mesh generation option multiplies the element number to 414,735 elements, almost 20 times more than "Extremely Coarse" option. However, the computational time needed to run the model using turbulence k-*ε* scheme is multiplied from 6760 s to almost 286,538 s, 40 times longer. The partial refinement mesh generation, used for models C, D, E, and F, is performed by selecting domains or regions of interest (ROI) to be discretized finer than other domains. Here, the observation chamber, T-pipe, and inlet pipe are selected as the ROI so that the velocity profiles and other important aspects in the region can be studied more accurately. This method aims to reduce the amount of mesh and computational time, and to avoid unnecessary data collection in non-ROI parts such as the diffusing chamber, drainage chamber, and converging chamber. The Partial Refinement method certainly opens numerous possible mesh arrangements. To limit the study, we focus on mesh refinement of the ROI only. The refinement of the observation chamber sub domain will enhance accuracy of the velocity profiles in the region needed for possible twodimensional modeling. On the other hand, the refinement of the inlet pipe and the T-pipe is expected to give accurate information on the amount of pump pressure needed in the

design. These three domains will be refined using "Extremely Fine" mesh or customized mesh, while the remaining body will be meshed using either "Normal" of "Finer" scheme. For Model C, the ROI are meshed using "Extremely Fine", while the remaining domains are meshed using "Normal". For Models D and E, the mesh for the test chamber domain of Model C is further refined by reducing the "maximum element size" from 0.0318 to 0.015 m. Model D is the only model that is computed using the HP desktop computer. It should be noted that the "Normal" size used in the non-ROI domains of models C, D, and E employs "maximum element size" that is about triple than that of the "Normal" size used in Model B4 with global mesh method. Consequently, although the ROI domains of Models C, D, and E, are finely meshed, the total number of elements of these models are less than that of Model B4. The difference in the "maximum element size" is due to the different "Calibration" option offered by COMSOL (COMSOL Inc., Burlington, MA, USA). The mesh configuration for Model F is the finest among models with partial refinement mesh. Here, the mesh configuration of the ROI used in Model D is maintained, but the remaining domains are meshed using "Finer" mesh. This refinement results in 355,730 elements, and the computation time is 145,860 s, only half of the time needed by model B4.

**Table 7.** This table lists CFD models used in this project and their associated number of elements and computation times. Only models A1 and A2 are executed using Laminar fluid model. All other models are executed using Turbulent model—RANS K-*ε* . All models are run in the HP ProBook, except Model D.


Shown in Figure 4 panel (b), results of Partial Refinement strategy where high mesh density ("Extremely fine mesh") is applied on the observation chamber, pipe inlet, and T-junction domains. The remaining domains are discretized using "Finer" mesh. The minimum and maximum element sizes employed in the observation chamber are 5 mm and 15 mm, respectively. The application would result in at least 10 computation nodes across the 15 cm wide observation chamber, which should be enough for accuracy.

The volumetric flowrate (*Q*) and mean velocity (*V*) in the observation chamber in the inlet pipe and observation chamber will be used to measure the accuracy of the solution. The volumetric flowrate across planes perpendicular to the long direction of the observation chamber must equal to that of the inlet. This flowrate can be obtained from the uniform inlet velocity *Vin* provided at the pipe inlet; *Q* = *Vin <sup>π</sup>* <sup>4</sup> *<sup>D</sup>*2, where *<sup>D</sup>* is the pipe diameter. The mean velocity in the observation chamber can be easily estimated using *V* = *<sup>Q</sup> Ach* , where the cross-section area of the chamber is known as *Ach* = 0.152 = 0.0225 m2. As the outlet at the base of the drainage chamber is assumed to have a pressure drop, the flowrate at this location should also be verified.

Second only to the mesh generation, the selection of appropriate physical modeling of the fluid is the most crucial aspect of CFD modeling. The selection between Laminar and Turbulent models is often determined by Reynolds numbers involved in the problem. The laminar threshold for fully developed flow in long pipe with circular cross sections is typically around 2300, while for an open channel flow is 500. Students have to understand that the change from laminar to turbulent flow is a gradual transition, and therefore it is

not practical to point out a single number to differentiate the two regimes. As pointed out by Lowe, the confusion over the critical Reynolds number is not new in the Fluid Mechanics community [49]. Nevertheless, smaller Reynolds numbers should be treated as turbulent when the complex geometry of the domains include circulation zones that are best captured using turbulent models. In the current water tank design, the maximum expected mean velocity occurring in the 15 × 15 cm square observation chamber is 60 cm per second, which results in the Reynolds number of 90,000, which puts the model under turbulent category. A long list of turbulent models is available in COMSOL Multiphysics software (COMSOL Inc., Burlington, MA, USA). For this project, we selected the standard k-*ε* model to serve our purpose. This suggested model is selected due to its popularity for industrial applications, easy convergence, and low memory requirement [50]. Alternatively, the analysis on such water flume can be performed using k-*ω* SST model [51]. However, a comparison study performed on a water flume with 5 m wide and 17 m long test section suggest that the outcomes from the two different Turbulent models do not show significant differences [52]. Heyrani et al. compare the performance of seven (7) turbulent models used in steady state modeling of a venturi flume. It was discovered that the k-*ε* model performs slightly better than the k-*ω* SST model [53].

Figure 5 shows the accuracy and computation times that are plotted against the mesh number for different models. The graph shows that Turbulent model with Partial Refinement (rectangular markers) produces the best accuracy among the three models (Laminar model with Partial Refinement, Turbulent model with Global Mesh, and Turbulent model with Partial Refinement). Moreover, its computing time is less than the Turbulent model with Global Mesh (circular markers). The computing time of the Laminar model with Partial Refinement (triangular marker) is the lowest, but the accuracy is also very low. The laminar model should only be used to make sure that the CFD model can be properly executed. High accuracy can be easily achieved by the Turbulent model with Partial Refinement, even with a low amount of mesh.

**Figure 5.** Convergence studies—data with dark labels indicate relationship between number of mesh and computation time, while the white markers indicate the relationship between the number of mesh and accuracy. The Laminar model shows the lowest accuracy and computing time. The Turbulent models can take 10 times longer time than the laminar flow, but its accuracy is very high.

Streamlines in the water tank when the mean inlet flow is 3 m/s are shown in Figure 6. The streamlines are obtained using "standard point controlled" method and using "Number of Points" entry of 80 points. The streamlines on the left and right are obtained when the computation modules used are Laminar model and Turbulent k-*ε* model, respectively. The streamlines from the Laminar model are parallel and straight in the observation chamber, but they are clearly absent from the circular flows in the drainage chamber, as well as in the diverging chambers that are demonstrated by the Turbulent model.

**Figure 6.** Three-dimensional streamlines obtained from an (**a**) Laminar model and (**b**) Turbulent model. The Turbulent model captures well the twin circulation zone occurring in the drainage tank and the collecting chamber. The Laminar model fails to capture the circulation zone despite of the use of large number of elements.

Observing the axial velocity profiles along the observation chamber, the Laminar model shows reasonable parabolic shapes (data are not shown here), however the amount of flow rate and mean velocity show deviation from the correct values. The profile of axial velocity along the observation chamber obtained from the Turbulent models are depicted in Figure 7. Each velocity profile is obtained at the same horizontal mid-section of the observation chamber, but at different distances "Y" from the starting of the chamber (taken as the interface between the converging chamber and observation chamber). The development of the velocity profile can be observed in this figure as the profile changes into a more parabolic shape as the Y is increasing. The velocity profiles demonstrate the typical characteristic of turbulent profile consisting of a turbulent core in the middle and laminar sublayer near the wall [45]. Note, however, that for the Model B4 (global mesh), the laminar sublayer is pronounced, while for Model F (partial refinement), it is very difficult to observe due to its small thickness. The accuracy of Model F in predicting the expected flowrate and mean velocity however is higher than that of Model B.

The velocity map and streamlines, viewed from the top of the flume, can be seen in Figure 8. The flow direction goes from top to bottom (diffusing chamber to the drainage chamber). Panel (a) on the left shows the velocity map of laminar model, while panels (b) and (c) are from the turbulent model, but with different mesh schemes. All models present reasonable straight parallel streamlines in the observation chamber. The Laminar model shows gradual increase in speed as it enters the observation chamber, as can be observed from the changing color. On the contrary, the Turbulent models both show consistent average velocity along the observation chamber. The Laminar model fails to produce the twin circulation zones that occur in the drainage chamber. These circulation zones are physically observed during real experimentation. The two Turbulent models clearly show these circulation zones with slight difference. The circulation zones by the Model B4 (Figure 8 panel (b)) are shown to be less chaotic than that of Model F (Figure 8 panel (c)).

**Figure 7.** Axial velocity profiles obtained at the horizontal mid-section of the observation chamber taken at different distances Y from the entrance obtained using (**a**) Global Mesh refinement and using "Normal" scheme (Model B4) and (**b**) Partial Refinement method (Model F).

**Figure 8.** Streamlines in the observation chamber and the diffusion chamber as well as in the collecting chamber produced by (**a**) Laminar model, (**b**) Turbulent model with Global Mesh "normal" refinement, and (**c**) Turbulent model with Partial Refinement method. The mean velocity of these models at the inlet is 3 m/s.

Vorticity is a local measure of the rotational motion of fluid particle relative to its own centroid, while the circulation can be considered as the measure of the rotational motion in global sense, relative to a distance reference point. After passing the observation chamber, water enters a short drainage chamber that is three times wider than the observation chamber but much deeper in the vertical sense. Here, due to the high speed in the observation chamber, the water shoots into the back wall of the chamber, but it also starts to fall due to the gravity. Hence, the stream splits into a pair of almost identical spiraling flows downward. Shown in Figure 9 are streamlines projected on a two-dimensional plane cut across approximately the middle height of the drainage chamber. On each panel, the streamlines demonstrate a pair of spiraling circles running in opposite directions. Note that the flow is coming from the observation chamber located, according to these pictures, relatively above the rectangular areas shown in the panels. The color map describes local

vorticity values, and the scales indicate that, on average, the two pairs have a similar amount but in opposite directions. In this figure, each panel represents a circulation zone in the chamber from the same Reynolds numbers, around 72 K (uniform inlet velocity of 3 m/s), but each is computed using different mesh density and element types used in the drainage chamber domain. On panel (a), the fluid model is laminar, and the elements used in the drainage chamber is tetrahedral. On panels (b), (c), and (d), the k-*ε* Turbulent model is used, but the number of elements used (only in the drainage chamber) are 5498, 14,681, and 61,052, respectively. On panel (b), only tetrahedral elements are used, but various elements are used for results shown in panels (c) and (d). It can be seen that the utilization of various elements results in smooth vorticity map and streamlines that are less discrete. The symmetry of the circular flow is captured when the mesh density is very high but, generally, panel (c) shows that essential flow structure is sufficiently captured. The maximum and minimum values on the scales indicate the amount of vorticity details that can be captured by the respected mesh density. On panels (b), (c), and (d), while the maximum and minimum vorticity are different, the average values inside and around the cores of the circulation area is around 4 to 5. On panel (c) and (d), zones with large vorticity up to 20/s can be captured, but these areas are concentrated near the walls. As is expected, the large vorticity occurs near the walls where large shear strain is expected to occur. On panel (a), the Laminar model results in a pair of seemingly symmetric circulation zones. The zones, however, appear to be shifted to the side walls. Moreover, the vorticity map shows some level of local spins, but the distribution is different from that shown in other panels. The laminar model certainly has failed to produce the desired outcomes, despite of the high mesh density that has been employed.

**Figure 9.** Panel (**a**) shows the circulation zone captured using laminar modeling; (**b**) Partial Refinement results in 5498 tetrahedron elements in the drainage tank, a minimum element quality of 0.6541, and an element volume ratio of 0.001275; (**c**) Global "Coarser Mesh" mesh with Turbulent model with 14,681 mixed elements, a minimum element quality 0.1307, and an element volume ratio 0.0021; and (**d**) Global Mesh refinement with "Normal" mesh with Turbulent model with 61,052 mixed elements in the drainage tank, a minimum element quality of 0.1251, and an element volume ratio of 9.8 <sup>×</sup> <sup>10</sup>−4.

#### **5. Fabrication and Manufacturing**

Welding was the most challenging process of the manufacturing activity in this project. All fabrication (cutting, folding, welding, assembly, etc.) had to be performed mostly

off-site at a private workshop due to the space limitation for hot works in our campus. Particularly, the welding of individual manifold was the primary challenge. The inlet and outlet manifolds needed to be designed, cut, welded, and assembled without the benefit of a fully equipped fabrication facility. Furthermore, the limitations of the available workshop dictated that tolerances could be held only down to 1/4 " (~6.35 mm) for most features of size and 1/2 " (~12.7 mm) for some larger features. The choice of 18 GA (0.049" or 1.27-mm) Austenitic Stainless Steel for the manifolds dictated that TIG (Tungsten static electrode Inert Gas shielded) welding was necessary to accomplish production. Alternative processes were ruled out for manifold production due to heat transfer concerns that would lead to excessive warpage of the fabricated sheet metal structures. Additional requirements of the design included the need for flush weld joints which would be less easily accomplished than other welding processes. Shown in Figure 10 one of the completed holding tanks (left) prior to being assembled with the inlet manifold (top right). The student who performed the welding used many cardboard templates (such as the one shown on panel (a) of Figure 10) to determine the cutting of the steel plate.

**Figure 10.** Panel (**a**) shows the cardboard model used in the making of the diffusion channel. Panel (**b**) shows the observation chamber being assembled, and panel (**c**) shows the final product of drainage chamber before being assembled into the main body.

MIG (Metal wire feed electrode Inert Gas shielded) welding was chosen for frame production. Mild Steel Square Tube with 1/16" wall thickness (~1.59 mm) was chosen for the frame and could withstand MIG welding with minimal concerns for warpage. The design did not require the welded joints to be flush with the welded materials. The steel square tubes that were welded together required a coat of rust protective paint due to the material properties that would cause rusting. A Miller Multi-Matic 220 multi-process welder (Miller Electric Manufacturing Co., Appleton, WI, USA) was used for all welding processes, selected due to its versatility of process selection and advanced user interface and arc control systems. The frame had undergone minor modifications from the original design to acclimate for the weight of the pump, transportation, and vibrations that would pass through the system. Jack screws, similar to stands, allowed the system to elevate off the wheels to limit vibration through the floor, level the entire system and to relieve stress from the wheels, as their specifications stated that they would only withstand 544.31 kg. The frame that the pump was set onto had its own set of jack screws that allowed it to be lifted off the frame for the rubber insulators to absorb vibrations produced by the pump. The 15 <sup>×</sup> <sup>15</sup> <sup>×</sup> 60 cm3 observation chamber was made of 0.635 cm thick transparent plexiglass. The corners are secured using both glue and bolts. To prevent leaking, rubber gasket is applied throughout. Finally, several holes are drilled on the upper edge of the chamber wall (not shown in the figure) to allow experiment platform to be secured onto.

#### **6. Experiment Results and Discussion**

The following experiments were conducted using the water tank after its fully built to test its capability:


#### *6.1. Velocity Measurement of Flow in the Observation Chamber*

The speed of the water flow can be controlled by adjusting the rotational speed or frequency of the pump's propeller using the programmable speed controller. In the first experiment, we wish to establish the relationship between the frequency of the propeller and the resulting water velocity in the observation chamber. The average velocity in the observation chamber is measured using a flow meter by Vernier (Vernier, Beaverton, OR, USA) (product number FLO-BTA) shown in Figure 11 panel (a). The sensor is designed for external flow measurement such as flow in open channels and rivers. The measurement range of this affordable sensor is up to 4 m/s, with a typical resolution of ±1.2 mm/s, and an accuracy of about 1 percent. A data acquisition card, LabQuest Mini, also by Vernier (product number LQ-MINI), is used to collect the data and transfer the data to a laptop through a USB cable. The dynamic data can be displayed by Graphical Analysis GUI software available for free by Vernier. The data can be exported to either Excel, Matlab, or Phyton, for further analysis and post processing. Typical data sampling in our experiment is 1 Hz.

**Figure 11.** The panel (**b**) displays the velocity reading by Vernier flow meter versus time when the pump operates at 40 and 50 Hz. The propeller is placed either near the floor of the chamber ("Bottom") or approximately in mid-height ("Middle"). Data suggests minimum difference by the two methods. The panel (**a**) shows the propeller being lowered in the "Middle" position in the observation chamber.

During one set of data collection, we seek to check the speed variability across the chamber's depth. Here, we perform the flow measurement by placing the sensor's propeller near (about 3 cm above) the chamber's floor and approximately at mid-height above the chamber's floor. The data collection for each trial was done at least for 40 s. On each trial, the frequency of the pump is increased from 10 Hz to 60 Hz using intervals of 10 Hz. This data collection results in insignificant variability of the time-averaged speed between the two locations. Examples of data collected over a period of 120 s for pump frequencies of 40 Hz and 50 Hz are shown in Figure 11 panel (b). The transient regime for about 25 s at the beginning of data collection can be easily identified in both data. The steady state regime maintains a constant average speed and indicates small fluctuations. Data with solid labels are taken for measurement at "Bottom", while data with open labels are from "Middle". Data indicate that the placement of samples in between these two locations would result in more-or-less identical outcomes.

#### *6.2. Tests on Flow Straightener Designs*

The second set of experiments involved velocity measurement in the presence of a flow straightener. The flow straightener functions to dissipate the turbulence water flowing into the observation chamber. The three (3) designs by students are shown in Figure 12; Models 1, 2, and 3 displayed in panels (a), (b), and (c), respectively. The cross section of all the models is 15 × 15 cm. Model 1 was constructed by PVC tubes of the same diameter (12.7 mm) cut at the same length (~10 cm) that are glued along its long sides. These cut pipes were then arranged in parallel, and bonded together on its long sides to form a channel. This model has the largest area of holes. Model 2 is a 3D-printed PLA plastic rectangular block of 15 <sup>×</sup> <sup>15</sup> <sup>×</sup> 5 cm3, furnished with 12-by-12 holes, each about 1.2 cm in diameter. Model 3 is a combination of parallel 0.5"-PVC (12.7 mm in diameter) pipes that are hold together using 3D-printed blocks with holes on their ends. Among the three models, Model 1 is the cheapest, but is very tedious to make, as each pipe has to be glued one by one. During the bonding process, the pipes can be placed in between two parallel walls. Model 2 is the most expensive, when the printing time is included in the calculation. While the SolidWorks (Dassault Systèmes SolidWorks, Inc., Waltham, MA, USA) CAD model for this model is easy to make, it took more than 12 h to print the model using a regular desktop 3D-printer. The thin walls separating the holes require the model to be printed with high density (high infills), and this prolongs the printing process. While involving 3D-printed blocks, Model 3 is inexpensive, as the holes on these blocks are well separated from one another, requiring less printing density compared to Model 2. Nevertheless, the arrangement of the pipes into these blocks is quite tedious. In terms of durability, it was quickly found that Model 1 disintegrated quite easily after several tests in the water. Model 2 also shows signs of plastic fibers detached from the bulk structure after several tests. Model 3 is the most durable, but its effectiveness is the lowest.

The data of flow velocity, taken when the pump's speed is 10 Hz, versus time for different models of flow straighteners, are presented in Figure 13. The blue, red, black, and green lines indicate water velocity without flow straightener, with Model 1, with Model 2, and with Model 3, respectively. On average, the velocity data for Model 1 (red line) is very similar to the case where no flow straightener is installed (blue line). Some degree of deviation from the flow with no straightener (blue line) is reflected by data using Model 2. Compared to other models, Model 3 shows the largest deviation from the flow with no flow straightener (blue line). Clearly, the Model 1 produces the least flow disturbance in terms of time-averaged velocity outcome. The expected outcomes can be attributed to the difference in the cross-section areas of the holes of each model. The largest area of holes for Model 1 is contributed by the holes by the pipes and the holes in between these pipes. Nevertheless, the Model 2 can be seen to produce a more stable velocity (black line) across the time compared to Model 1 and Model 3.

**Figure 12.** Models of flow straightener used in the flow experimentation—(**a**) Model 1 is made of cut PVC pipes that are glued together in parallel arrangement. (**b**) Model 2 is a full 3D-printed apparatus designed in SolidWorks. (**c**) Model 3 is an assembly of PVC pipes fitted into 3D-printed block with holes.

**Figure 13.** This figure shows mean velocities from the flowmeter for various flow straightener models taken for a period of 120 s when the pump speed is 10 Hz. The blue line is the velocity of water flow without flow straightener. The velocity after the installation of Model 2 shows stability compared to other models.

Time-averaged (over the steady state period) water flow velocity obtained from the flow meter for various pump's frequencies are presented in Figure 13. The data points are labeled using white circles (no flow straightener), white square (Model 1 of flow straightener), and white triangle (Model 2 of flow straightener). A linear relationship of

$$\nabla = 0.014 \times Fr - 0.042 \left( R^2 = 0.9999 \right) \tag{4}$$

was discovered for data obtained without the flow straightener. Here, *V* (m/s) and *Fr* (*Hz*) are the mean velocity and pump's frequency, respectively. The non-zero intercept suggests the minimum frequency (associated with a minimum water velocity) needed to turn the propeller of the flow meter of about 3 Hz. The known linear relation between the pump's frequency and average velocity establishes the formulation of *V vs*. *F*. for future experimentation. Velocity data obtained using the flow straightener Model 1 shows similarity to the data obtained without the flow straightener. On the other hand, Model 2 produces velocity data that deviates at high pump's frequency. The linear relationship is also quantified as *<sup>V</sup>* <sup>=</sup> 0.0166 <sup>×</sup> *Fr* <sup>−</sup> 0.0663 *R*<sup>2</sup> = 0.9951 .

On a separate, simple experiment, two students measured the water velocity by means of a Lagrangian approach. Here, the time needed by a ~1 <sup>×</sup> <sup>1</sup> <sup>×</sup> 2 cm<sup>3</sup> cork to travel for 50 cm is measured using a manual stopwatch. This measurement results in smaller velocity compared to ones measured using the propeller. These Lagrangian velocities are displayed in Figure 14 using dark cross legends. The lower estimation can be attributed to the drag produced by the weight of the cork, as well as accuracy of the manual data collection. Nevertheless, the manual measurement using cork provides a simple and quick validation of the Eulerian method employed in the automatic measurement using digital flow meter.

**Figure 14.** This figure shows the comparison of time-averaged velocity data versus pump's frequency. The data are obtained using Models 1 and 2 of the flow straighteners. The data from measuring the water velocity using floating cork are presented using cross symbols. Data from Model 1 shows minimal deviation from the reference data. The cork data underestimates the reference data, but it gives good estimate for a quick verification of the flow meter.

#### *6.3. Flow Visualization Using Dye*

The effects of placing the flow straightener in the water stream was examined by means of dye injection. While this method is simple and affordable, students performing the tasks can learn valuable lessons. For example, students realize that the fluid density had to be carefully chosen, as a fluid that was too light would rise to the top (float) and the fluid mixture was too dense then it would sink to the bottom. The fluid was also mixed with food coloring to make it easier to observe. After running several tests using hand soap, oil, sugar, and honey, it was concluded that honey diluted with a small amount of water gave the best results. These tests were run by using a syringe with the mixed fluid. Figure 15 shows the difference in flow behavior through the utilization of injected dye. The left picture (panel (a)) shows the flow turbulence when there is no flow straightener inserted in the chamber. Panels (b) and (c) show the resulting flow behavior when Model 2 is used. Here, the flow is almost straight and uniform with minimum turbulence characteristic.

**Figure 15.** The flow behavior at about 10 cm/s. was observed using color dye injected by a syringe. The left panel (**a**) shows the dispersion of the dye when no flow straightener was used. Turbulence characteristic can be seen from the dyes. Panels (**b**,**c**) shows the flow characteristic when Model 2 is used as flow straightener. Straight streamline is easily observed in this case.

#### *6.4. Measurements of Drags*

Objects exposed to fluid flow experienced both drag and lift forces. Being able to demonstrate this phenomenon for objects in the water flow is one important objective of the water flume project. Experimentation on drag reported below was performed by students in Thermal Fluid Lab course on Fall 2019. The drag force acting on objects, such as long prismatic cylinders and bars, is registered using a force sensor from Vernier (Go Direct ® Force Sensor—GDX-FOR) that uses strain gauge to measure as small as ±0.1 N up to ±50 N of axial force, sufficient range for our fluid application. A sampling rate of 50 Hz, out of a possible maximum of 1000 Hz, is used in our application. This device can be either directly connected to a laptop using a USB cable, or wirelessly using a Bluetooth to a mobile device. The stream data can be displayed on a laptop and recorded using a software Graphical Analysis available for free from Vernier. The force sensor can also be connected to the data acquisition card by Vernier—LabQuest Mini. This allows simultaneous measurement of multiple sensors with other data such as water speed, temperature, and lift. The ~7 <sup>×</sup> <sup>4</sup> <sup>×</sup> 5 cm3 sensor device is placed on ~15 <sup>×</sup> <sup>15</sup> <sup>×</sup> 1.5 cm3 (1.36 kg) aluminum sled supported by four smooth wheels (rollers) that rest on the walls of the observation chamber (Figure 16). The sled was machined to include a hollowed pin on its top to secure the force sensor. A hole was drilled at the center of the sled to allow samples to be secured using bolt and nut. The other end of the bar is free. The force registration is performed by aligning the sensor's hook to the flow direction and securing it to a reference point fixed to the chamber's wall. The four smooth wheels reduce the contribution friction caused by the sled's weight and maximize the registration of the hydrodynamic drag by the sensor.

**Figure 16.** The figure on the left panel (**a**) shows the force sensor supported by the aluminum platform (sled) with its wheels placed along the wall of the observation chamber. The panel (**b**) on the right shows the complete setup of the drag sensor secured on the observation chamber with the flow meter placed upstream to allow simultaneous measurements of flow speed and drag forces. The water flows from left to right.

The drag measurement has been conducted for several objects: plastic cylinders, 3Dprinted cylinders, solid aluminum cylinders, and square bars, as well as rotating bladeless

turbine models. A preliminary attempt using a 3D-printed fixture to hold hollow plastic cylinder samples and 3D-printed samples resulted in promising data. However, it also indicated the need for rigid and sturdy support and samples that do not float. The force data is converted into Drag Coefficient, *CD* using the following formulation [43,44]:

$$\mathcal{C}\_D = \frac{2\mathcal{F}\_D}{\rho \nabla^2 A} \tag{5}$$

where *FD*, *ρ*, *V*, and *A* are the time-averaged drag force, water density, time-averaged water velocity, and frontal area of the cylinder sample, respectively. The frontal area only accounts for surface perpendicular to the water flow. Only data from the 16 mm diameter aluminum cylinder and square aluminum bar are presented in this report.

Shown in Figure 17 is an example of raw data of the force obtained from the force sensor for a period of 120 s when the data is taken at 50 Hz sampling frequency. The average velocity of the water flow here is about 59.4 cm per second. The oscillation of drag force is expected to be caused by the unsteadiness of the flow and the vortex shedding that occurs when fluid flow passes the object. The average force is about 0.38 N, and its standard deviation is ±0.072 N.

**Figure 17.** This graph shows typical force data collected during a test. For this data, the time-averaged velocity is 0.594 ± 0.01 m/s (not shown), while the time-averaged force is 0.386 ± 0.072 N. The sample is 16 mm cylinder bar using data sampling of 50 Hz.

The drag force data for each sample was collected for four different pump speeds: 0.35, 0.50, 0.60, and 0.66 m/s. All measurement was obtained when steady state had been reached. Figure 18 shows monotonic and non-linear increase in the drag forces for the samples in the range of the applied flow speed. The two samples demonstrate a very similar trend. The error bars show the standard deviation over the data population of speed and force measurement. The error percentage presented by the flow meter is less than that by the force sensors. A consistent increase in the deviation from the mean of the force data is observed when the flow speed increases. This is expected, as the force should be proportional to the square of flow speed. Data from the flow meter indicate consistent deviation throughout the range of the flow speed.

**Figure 18.** The mean velocity vs. mean force shows monotonic increase in the drag force with the increase in applied velocity. Error bars show standard deviations that increase with the velocity. The errors from the flow meter are relatively smaller compared to that from the force sensor.

Figure 19 shows the plot of the average drag coefficients for increasing Reynolds numbers. The average *CD* here is calculated using the time-averaged velocity and timeaveraged force taken at a specific pump's frequency. The error bars indicate the "maximum" and "minimum" values of the drag coefficients, calculated using the maximum and minimum velocity and force data reading. The large fluctuation of force data at high velocity is reflected here as large range of "maximum" and "minimum" coefficients. Nevertheless, the average drag coefficients for both cases remain relatively constant in the presented range of Reynolds numbers. The average drag coefficients, over the given Reynolds numbers, are 1.48 ± 0.25 and 1.58 ± 0.13 for the cylinder and square bar, respectively. The plot of *CD* for cylindrical bar as a function of Reynolds numbers for 0 <sup>≤</sup> *Re* <sup>≤</sup> 104, proposed by Tang et al. [54], is also shown in this figure. The expression is stated here for clarity:

$$C\_D = 10^{\left(0.0041 \times \left(\log\left(Re\right)\right)^3 + 0.0853 \times \left(\log\left(Re\right)\right)^2 - 0.65 \times \log\left(Re\right) + 1.05\right)}\tag{6}$$

The function is obtained by regression analysis of *CD* data on cylinder bars obtained from various data collection. The dash line in Figure 19 shows the constant *CD* number for a long cylinder that is typically published in textbooks. The presented comparison shows that the *CD* values obtained in our lab can be considered reasonable. For *CD* of the cylinder bar, the deviation from both published *CD*s is large when the Reynolds number is low, but the error is minimized as the Reynolds number increases. The outcome certainly warrants further refinement of the experiment apparatus and more data collection at various Reynolds numbers.

**Figure 19.** Average drag coefficients (*CD*) for the circular cylinder and square bar for various flow speeds (Reynolds numbers). The error bars indicate the "maximum" and "minimum" *CD* values calculated using the extreme values of speeds and forces obtained during the sampling period of 120 s. The average *CD* values across the range of Reynolds number for the cylinder and square bar are 1.48 and 1.58, respectively. The data marked by dark circles are obtained using formulation proposed by Tang et al. [54], and the dash line represents the constant known *CD* typically presented in textbooks.

Alternatively, the drag coefficients can be obtained by first plotting the drag force data against the square of the time-averaged speed, as shown in Figure 20. Using standard linear regression analysis, the linear relationship between *FD* and *<sup>V</sup>*<sup>2</sup> can be obtained. For the presented samples, the following relations are obtained *FD* <sup>=</sup> 1.1623 *<sup>V</sup>*<sup>2</sup> and *FD* <sup>=</sup> 1.3449 *<sup>V</sup>*<sup>2</sup> for the cylinder bar and square bar samples, respectively. The high *<sup>R</sup>*<sup>2</sup> values for the two samples indicate a strong fit in the regression process. It also indicates that the linear relation is appropriate, as the *CD* is expected to be constant in the applied range of Reynolds numbers. Note that the zero intercept in this case is not automatic, and it must be enforced as the zero force should be correlated with zero velocity. Eventually, the drag coefficients can be calculated from the gradients of the lines using the following formulation, which can be derived from (4)

$$
\hat{\mathbf{C}}\_D = \frac{2 \text{ K}}{\rho \text{ }^\circ \text{A}} \tag{7}
$$

where *<sup>K</sup>* is the gradient of the linear relationship between *FD* and *<sup>V</sup>*<sup>2</sup> . The coefficients found using this method are very similar to the average coefficients, however this method results in coefficients that are closer to the known values from textbooks [43–45].

**Figure 20.** Linear correlation between the drag force and *<sup>V</sup>*<sup>2</sup> (square of the velocity) can be used to obtain the drag coefficients. The high *R*<sup>2</sup> numbers suggest that the linear regression fits the data very well. The gradient for the square bar data (empty square markers) is slightly higher than that of the cylindrical bar (circle markers).

Presented in Table 8 is a list of the *CD* values obtained in our experiments compared to those presented in textbooks. The 2D and 3D refer to the two-dimensional and threedimensional assumptions of the immersed bodies. The *CD* refers to the averaged coefficient of drag. This *CD* is obtained by simply taking the average of all *CD*s calculated using different mean velocities and their corresponding drag forces. The alternative *C*ˆ *<sup>D</sup>* is obtained by first plotting the drag force versus the square of the mean velocities (*FD vs*. *<sup>V</sup>*<sup>2</sup> ) and then computing the gradients of their linear relationship. The last column of the table shows the *CD* value when the cylinder is considered as a three-dimensional (3D) body. The length-to-diameter ratio of our sample is *<sup>L</sup> <sup>D</sup>* ∼ 7, and the corresponding *CD* value is in between 0.74 and 0.82 [44]. Compared to the 2D version of the known values, both methods result in drag coefficients within reasonable margin errors for the given experiment setup and range of applied Reynolds numbers (5000 <sup>≤</sup> *Re* <sup>≤</sup> 104). The known 2D *CD* is stated to be valid for 103 <sup>≤</sup> *Re* <sup>≤</sup> 104 (laminar regime), while the 3D *CD* is valid for *Re* <sup>≥</sup> 104 (turbulent regime) [44].

**Table 8.** List of the *CD* values obtained in our experiments compared to that presented in textbooks. The 2D and 3D refer to the two-dimensional and three-dimensional assumptions of the immersed bodies. The *CD* refers to the averaged coefficient of drag. The *C*ˆ *<sup>D</sup>* is obtained using the linear regression analysis of *FD vs*. *<sup>V</sup>*<sup>2</sup> .


Given the limitation of the support system and the accuracy of the sensors, the drag coefficients obtained in this experiment can be considered reasonable. The close walls presented by the observation chamber introduce effects on the vortex trail, hence the drag forces. Numerical studies at low Reynolds numbers showed that the wall increases the drag coefficient [55,56]. Experimentation at high Reynolds numbers (above 3000) on a cylinder, flat plates, and a square bar indicates the effects of walls on the pattern of the vortex-shedding [57]. The presence of the flow meter upstream of the sample, as mentioned previously, should have affected the force reading. Simultaneous flow rate reading is not necessary, as the relationship between the pump's frequency and flow speed has been known. Alternative supporting systems for the bar, for example pin support or fixed support on both ends, should be studied in future. Lastly, the absence of flow straightener needs to be considered, as it is known to affect the flow characteristics in the observation chamber.

#### **7. Budget**

Table 9 shows major items purchased for the flow tank, totaling approximately USD 5550. This is slightly higher than the USD 3500 cost of a custom-made water tunnel with a 30 cm wide and 180 cm long test section reported by Darfler and Tsai [2]. Our budget matches well with a flow loop with similar scale testing chamber reported by Northern Arizona State University [58]. However, this system can only perform laminar flow due to the small pump size. The centrifugal pump, as is expected, is the most expensive item in our project. However, its cost is still less than 22% of the total budget. The cost for the speed controller is also modest. The remaining cost covers the materials, such as stainless steel plates for the tank body, an acrylic sheet for the observation chamber, various pipes and valves for connectors, and steel tubes for the supporting frames. Gas supply for the welding takes less than USD 1000. The total cost listed here does not include the labor cost for the welding, assembly, design, transportation, wiring, and painting. All price is based on 2018 market values. An up-to-date pump price may be found on the manufacturer's website [59].

**Table 9.** List of major components of the flume, amount and its 2018 price. Note that the labor cost is not included here.


#### **8. Conclusions**

The design and manufacturing process of an educational water flume has been shown to provide a learning platform for engineering students. The project is suitable for a senior design capstone project that would take two (2) semesters to complete. The design process has provided complex hands-on exercise materials for students desiring to master engineering fundamentals and skills such as computer aided designs, computer analysis of structures and fluid flows, electrical circuits, and experimentation. The open flume can be an expensive investment, but its versatility truly provides avenues for multi-disciplinary creative research. For a small-scale water flume, the required pump is governed by the height of the observation chamber and desired flow rate for the experiments. Welding of the tank body and supporting frames take the most challenging step and must be facilitated well. Professional welding and manufacturing, when budget is allowed, should be preferred to obtain a quality product, as rusting can become imminent.

**Funding:** The construction of the flume was funded by the Vibration Institute through its Academic Grant program for 2018. The Capstone Senior Design Project was partially funded by the Division of Engineering, State University of New York (SUNY) at New Paltz, New York.

**Institutional Review Board Statement:** Not Applicable.

**Informed Consent Statement:** Not Applicable.

**Data Availability Statement:** All essential data have been presented in the paper. Additional data are available on request from the corresponding author.

**Acknowledgments:** The authors wish to thank many alumni of the Division of Engineering who have contributed to this project: Matthew Bogue, Peter Demertzis, Jaime De La Vega, Roneil Harris, Apurba Sharma, Seth Pearl, Jessy Li, Michelle Wong, Kieran Cavanagh, Kira A. Solano, Christopher M. Rodi, Hadi Rabadi, Maritsa Rehma, and Istvan "Steve" Benyei. Many thanks for the continuing support from faculty colleagues during the project: Yi-Chung Chen, Faramarz Vaziri, Eric Myers, and Ashwin Vaidya as well as the technical lab support from Anthony Denizard.

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

#### **References**


## *Review* **An Introduction of Droplet Impact Dynamics to Engineering Students**

#### **Sara Moghtadernejad 1,\*, Christian Lee <sup>1</sup> and Mehdi Jadidi <sup>2</sup>**


Received: 1 June 2020; Accepted: 1 July 2020; Published: 2 July 2020

**Abstract:** An intensive training course has been developed and implemented at the California State University Long Beach based on 8 years of experience in the multiphase flow area with the specific focus on droplet–solid interactions. Due to the rapid development of droplet-based equipment and industrial techniques, numerous industries are concerned with understanding the behavior of droplet dynamics and the characteristics that govern them. The presence and ensuing characteristics of the droplet regimes (spreading, receding, rebounding, and splashing) are heavily dependent on droplet and surface conditions. The effect of surface temperature, surface wettability, impact velocity, droplet shape and volume on droplet impact dynamics, and heat transfer are discussed in this training paper. Droplet impacts on moving solid surfaces and the effects of normal and tangential velocities on droplet dynamics are other topics that are discussed here. Despite the vast amount of studies into the dynamics of droplet impact, there is still much more to be investigated as research has expanded into a myriad of different conditions. However, the current paper is intended as a practical training document and a source of basic information, therefore, the scope is kept sufficiently broad to be of interest to readers from different engineering disciplines.

**Keywords:** droplet impact; undergraduate education; applications of fluids

#### **1. Introduction**

This paper summarizes an intensive course on the subject of multiphase flow with the focus on both the theoretical and practical knowledge on droplet–solid interactions. More specifically, the goal of this training manuscript is to provide students with in-depth practical knowledge of droplet dynamics due to its numerous industrial and scientific applications. The dynamics of droplet impacts on solid surfaces has been studied for years starting with the work of Worthington [1,2]. To this day, droplet impact dynamics are widely studied as further improvements in understanding of droplet behavior due to the variety of their industrial applications. In particular, droplet impacts are imperative to industrial processes employing spray coating and painting, spray cooling, inkjet printing, combustion engines, and anti-icing characteristics of critical industrial components such as aircraft surface, powder lines, and wind turbines [3–7]. A brief discussion for some of these usages is presented below.

One of the prominent applications of droplet studies is aircraft icing, which refers to creation of ice on the surface of flying objects. This usually happens at the presence of super-cooled water droplets and below icing temperature. Those droplets are formed as droplets with various size merge in clouds, or when the falling snow melts as it passes through a warmer layer (weather inversion). In either case, they are very unstable, and any disturbance will cause ice formation. During a flight, a fraction of super-cooled droplets tracks through the airflow and impact the aircraft surface. In most cases, the low energy droplets freeze upon impact, whereas the ones with higher energy flow along the surface until

their energy is depleted and then freeze. Consequently, the resulting ice changes the surface structure of the aircraft and this reduction in the aerodynamic of the wing increases the fuel consumption [8–13]. As such, a detailed understanding of this process is much needed to engineer more efficient wings.

Other application of droplets is in direct injection combustion engines, where the behavior of fuel droplet impingement on the piston and surrounding cylinder walls is heavily studied [14–18]. The formation of the wall-fuel film and the improvement in the efficacy of combustion engine promote a wide variety of research into droplet dynamics [14–18]. The properties of superhydrophobic surfaces are of major interest in aerospace and power industries for their anti-icing characteristics, to prevent ice formation on wind turbine blades, power-lines, or aircraft wings [11,19–21]. Methods improving the repulsion of water droplets off surfaces have been heavily studied. Such methods involve variable surface inclination and application of superhydrophobic coatings [4,22–26].

Another interesting application of droplet dynamics is in the thermal spray process, where molten or semi-molten metal and ceramic particles are deposited on a substrate to generate various types of coatings such as thermal barrier, wear resistance, and corrosion resistance ones [27–32]. During the thermal spray process first, a heated gas is created via chemical combustion or electrical energy in a torch. Then the gas is used to melt the coating powder or wire into droplets and accelerate them toward a substrate, where the particles generate a splat. Finally, a coating is formed via multiple layers of splats. In that regards, most research studies the effect of droplet's inflight behaviors such as velocity, temperature, and trajectory, and their interaction with the substrate's temperature and roughness on the process of droplet solidification and coating structure [27–31,33–38].

Yet another importance of fundamental study in the droplet impact field is its application in spray cooling, where an array of small droplets is applied to a heated surface as a cooling mechanism to enhance its heat transfer [39–44]. In this operation, the cooling effectiveness is strongly influenced by fluid properties and droplet's size and velocity. Amongst the numerous usages here are dermatological operation, fire protection, and cooling of hot surfaces like hot strip mill and high-performance electronic devices. For example, cryogenic spray cooling is selectively directed to pre-cool human skin in laser treatments and hair removal procedures. Similarly, in steel strip casting, a jet gas with water droplets is guided to cool a high temperature (up to 1800 K) steel surface that shapes the final microstructure optimization. All in all, new developments in spray cooling technology demand improving the heat transfer rate, while maintaining uniform heat removal, and preventing temperature overshoot. They also require uniform operating temperatures maintenance, the removal of high heat flux, or the adaptability to changes in heat flux [39,42,45].

Furthermore, droplet dynamics are widely used in ink-jet printing as it involves the generation and deposition of small droplets, usually containing colorants, onto a substrate in certain patterns. Main challenges in the process are when droplets bounce back or spread unevenly on the paper. The development in understanding droplet dynamics is thus motivated by the vast applications of ink-jet printing, where improvements in droplet positioning, volume, and directionality are the main concerns [46–48]. Note that, the input substance (ink) can be a tiny liquid of "smart" material. As such today's applications of ink-jet technologies include printing solar cells, medical sensors, and electronic circuits. Particularly, ink-jet printing techniques have rapid developments in medical industry, where the manipulation of small amounts of liquid is essential in cell culture growth and pharmaceutical drug production [49–51].

Lastly, there has been rapid developments in microfluidic technologies, where liquid drops are handled on nano or micro scales. Microfluidic technologies have heavily impacted the biomedical, environmental, food, and chemical industries [52–54]. Other notable industrial applications of droplet impacts include quenching of aluminum alloys and steel, fire suppression, incinerators, soil erosion, and crop spraying [3,55–58].

#### *Key Dimensionless Numbers in Droplet Impact Dynamics*

To facilitate the characterization of droplet dynamics several dimensionless numbers have been introduced in the literature. The main dimensionless numbers are given in Table 1 [3,59,60], where ρ*<sup>l</sup>* denotes the liquid density, *d*<sup>0</sup> is the droplet diameter, *u*<sup>0</sup> is the droplet impact velocity, σ is the surface tension, μ*<sup>l</sup>* is the liquid viscosity, *g* is the gravity, and *t* is time. Furthermore, in this table, *dcl*, *cl*, *kl*, *T*0, *hlv*, *Q*, and *Mevap* denote the diameter of the wetted region, liquid specific heat, liquid thermal conductivity, initial temperature, enthalpy of vaporization, total heat that is transferred to the droplet, and the evaporated mass at a given time, respectively. Subscripts *l* and *s* also stand for liquid and solid, respectively. It is worth mentioning that Table 1 only covers the main dimensionless numbers, while additional dimensionless numbers can be defined for the contact angle of the droplet on the surface, the boundary layer thickness, vapor properties, etc.


**Table 1.** Key dimensionless numbers to study droplet impact dynamics.

#### **2. Droplet Impact on Static Solid Surfaces**

Experimental investigations suggest the presence of six possible regimes for droplet impact on dry surfaces, including deposition, receding breakup, rebound, and splashing (see Figure 1). The characteristics of these regimes after droplet impact has been deeply investigated [3,57,59,61–64]. The presence, or lack of each regime is dependent upon multiple conditions of the impacting droplet and characteristics of the surface. Such properties include impact velocity, droplet size, liquid surface tension and viscosity, surface temperature, wettability [3,4], etc.

**Figure 1.** Droplet impact regimes on dry surfaces (reproduced with permission from [65]).

The expansion of the droplet after its impact on the surface is known as the spreading regime. This phenomenon is observed at relatively low Weber number values and low surface temperatures [3,63,66]. One of the most important parameters for characterizing the spreading regime during the droplet impact is the maximum spreading diameter of the droplet. The maximum spreading diameter of the droplet is imperative when heat transfer is a concern as the maximum spreading defines how much surface area is available for heat transfer [66]. Several characteristics of both the droplet and surface can affect the maximum spreading, such as droplet impact velocity and surface wettability [67–70]. Initial droplet velocity heavily impacts its spreading. As the droplet impacts on the surface, the kinetic energy in the vertical direction is transferred to the radial direction, promoting the spreading of the droplet. Lann et al. [69] described in their work that droplets spreading behavior is overseen by the conversion of kinetic energy into surface energy or dissipated heat. Different models have been developed to predict the maximum spreading ratio, *Smax* = *dcl*,*max*/*d*0, as a function of the impact parameters. Some of the most commonly used formulas in the literature are as follows:

Jones [71]

$$S\_{\text{max}} = \sqrt{\frac{4}{3} R e^{1/4}}\tag{1}$$

Chandra and Avedisian [72]

$$\frac{3\mathcal{W}\epsilon}{2Re}S\_{\text{max}}^4 + (1 - \cos\theta)S\_{\text{max}}^2 - \left(\frac{\mathcal{W}\epsilon}{3} + 4\right) = 0\tag{2}$$

Asai et al. [73]

$$S\_{\text{max}} = 1 + 0.48 \mathcal{W} \varepsilon^{0.5} \times \exp\left[-1.48 \; \mathcal{W} \varepsilon^{0.22} \; R \varepsilon^{-0.21}\right] \tag{3}$$

Pasandideh-Fard et al. [74]

$$S\_{\text{max}} = \sqrt{\frac{We + 12}{3(1 - \cos \theta\_A) + 4\left(We / \sqrt{Re}\right)}}\tag{4}$$

Mao et al. [75]

$$
\left[\frac{1}{4}(1-\cos\theta)+0.2\frac{\mathcal{W}e^{0.83}}{\mathcal{R}e^{0.33}}\right]\mathcal{S}\_{\text{max}}^3 - \left(\frac{\mathcal{W}e}{12}+1\right)\mathcal{S}\_{\text{max}} + \frac{2}{3} = 0\tag{5}
$$

Ukiwe and Kwok [76]

$$\left(\left(\mathcal{W}e+12\right)\mathcal{S}\_{\text{max}} = 8 + \mathcal{S}\_{\text{max}}^3 \times \left[3\left(1-\cos\theta\_A\right) + 4\frac{\mathcal{W}e}{\sqrt{\mathcal{R}e}}\right] \tag{6}$$

where θ*<sup>A</sup>* stands for the advancing contact angle.

The receding regime involves the contraction of the droplet once it has reached the maximum spreading. Droplet receding is due to the liquid surface tension, which pulls the droplet together [3,66]. The receding time depends on several factors including the surface shape and the liquid surface tension [67,77]. The suppression of this regime has also been observed under various conditions. For example during the impact of liquid fuels such as decane, ethanol, and tetradecane due to their low surface tension, the contraction phase is suppressed [77]. This regime has also been shown to be heavily affected by surface temperature. Receding can be hindered or even vanished by sufficiently low surface temperatures. For instance, during the water droplet impact on ice films, once the droplet reaches maximum spreading, the contact interface between the droplet and the ice film freezes, preventing the formation of a receding regime [78–82].

After receding, the droplet may show rebounding depending on various parameters such as surface wettability, temperature, and kinetic energy. In the rebound regime, if the droplet has a high contact angle and sufficient kinetic energy, it may bounce off the surface. The rebound regime has also been shown to be absent under various conditions. For instance, upon impact on a super-cooled superhydrophobic surface, water droplet freezes and adheres to the surface as it spreads, therefore, no rebounding occurs [83].

The appearance of the splashing regime is also dependent on various conditions of the droplet and surface. During impact with high enough energy, the droplet hits the surface and disintegrates into secondary droplets. Splashing can be induced by changing the orientation of the surface, such as having spherical or inclined surfaces, or impact on a moving liquid film [84–86]. There are several types of splash such as corona, prompt, and fingering within the splashing regime [77,87–89]. In corona splashes (see Figure 1) the outer rim of the lamella lifts off the surface forming a crown. From the crown, the lamella breaks apart forming secondary droplets. In prompt splashes, the droplet disintegrates, and secondary droplets form immediately after impact. Finally, in fingering splashes, protrusions extend from the droplet, eventually disintegrating and forming secondary droplets [77,87–89]. Various boiling regimes have also been experimentally studied at super-heated temperatures. This causes different splash phenomena, such as boiling-induced breakup or boiling-induced splashing, due to nucleate boiling in the contact interface [66].

As mentioned earlier in this manuscript, the regimes and phenomena observed during droplet impact significantly depend on a variety of droplet/surface conditions. Here we aim to describe the key conditions affecting the impact dynamics of droplet on solid substrates. As such the article is divided into two sections describing the droplet impacts on (1) stationary solid surfaces and (2) moving solid surfaces. Noting that due to the vast number of characteristics defining the behavior of droplets only the major ones will be discussed in this article; namely impact velocity, droplet shape and size, surface wettability, and temperature.

#### *2.1. E*ff*ect of Surface Temperature*

Surface temperature is an important characteristic affecting the dynamics of impacting droplets, particularly when heat transfer is a concern. In the case of a hot, dry solid surface, the droplet impact outcomes are classified into several regimes: evaporation, nucleate boiling, foaming, transitional boiling, and film boiling. In the nucleate boiling regime, which occurs at relatively high surface temperatures, the droplet is in direct contact with the surface and vapor bubbles are formed at various isolated nucleation sites. These bubbles rise and the droplet will eventually boil off (see Figure 2b) [67]. In the foaming regime, which is a subcategory of nucleate boiling, the entire drop starts to foam. The vapor bubbles grow much larger in this regime while no separation from the liquid–gas interface

and no coalescence are detected. In the transition boiling, due to high wall temperature, the generation rate of the vapor bubble increases quickly. Owing to bubbles coalescence, a vapor layer is formed over some portions of the area between the drop and the surface, while the rest of the drop wets the surface. In this regime, liquid layers frequently collapse, therefore, this regime is very unstable and secondary droplets are also generated [90]. Further increasing the temperature changes the regime to film boiling. In this regime, a vapor layer forms preventing the complete contact between the liquid and the surface [67]. In the film boiling regime, the droplet levitates on a vapor layer, as shown in Figure 3. This phenomenon is heavily studied and called the Leidenfrost effect [91]. The Leidenfrost effect also promotes rebounding of the droplet without disintegration into secondary droplets.

**Figure 2.** Outcomes of droplet collision with a hot, dry solid surface: (**a**) evaporation, nucleate (**b**) boiling, (**c**) foaming, (**d**) transition, and (**e**) film boiling (reproduced with permission from [90]).

**Figure 3.** The Leidenfrost phenomenon forms vapor layer preventing droplet from contacting the surface.

One important characteristic of heat transfer in droplet dynamics is the evaporation time. When the surface temperature is within the nucleate boiling regime, it was observed that as surface temperature increases, the evaporation time decreases due to an increase heat transfer. During the transitional regime, the evaporation time is variable due to uneven contact between the liquid and surface. However, in the film boiling regime, the evaporation time sharply increases due to the Leidenfrost effect. When the droplet enters the film boiling regime, it is no longer in contact with the surface as the vapor layer separates the two. This vapor layer acts as an insulating layer, slowing down evaporation [67]. Further increasing the temperature past the Leidenfrost effect will decrease the evaporation time [92].

Liquid surface tension has significant influence on the dynamic Leidenfrost temperature. The dynamic Leidenfrost temperature is the minimum surface temperature at which the impacting droplet bounces without splashing, in other words the minimum temperature to induce the Leidenfrost effect [93]. Chen et al. [93] studied the effects of surfactants on the dynamic Leidenfrost temperature. In these sets of experiments alcohol surfactants, octanol and ethyl-hexanol, were added to water droplets causing a reduction in surface tension of the droplet. As a result, the maximum spreading diameter increased during the spreading regime. Additionally, it was observed that the time to reach maximum spreading was decreased. Chen et al. [93] found that by adding surfactants the dynamic Leidenfrost temperature increased. The addition of surfactants promotes splashing and prevents the formation of vapor layer. The reduced surface tension promotes a thinner lamella during spreading, this makes it easier for vapor bubbles to burst out from the boiling film and break the vapor layer, preventing the Leidenfrost state [93]. The reduced surface tension also makes it difficult for

vapor bubbles to coalesce and form a vapor layer as a result of the lower bubble departure diameter. The dynamic Leidenfrost temperature can also be raised higher by increasing the concentration of surfactants in the droplet due to increased reduction in surface tension. Surface wettability is another parameter that affects the Leidenfrost point in a way that more hydrophilic surface leads to a higher Leidenfrost temperature [67,90].

Cooled and super-cooled surface temperatures also have an intense effect on droplet dynamics. At sufficiently cool temperatures, droplets begin to nucleate and freeze. The sessile droplet undergoes a complex solidification process, typically split into five stages [83,94,95], which is similar to the freezing of suspended water droplets. First, the droplet is cooled from its initial temperature to temperatures below the equilibrium freezing temperature during a cooling stage. Second, the droplet experiences a nucleation stage, where ice crystal nucleation occurs. Third, rapid crystal forms from the nucleation points driven by supercooling during the recalescence stage until it reaches the equilibrium temperature. Fourth, crystal growth is driven by heat transfer until the droplet is completely frozen during the freezing stage. Lastly, during the solid cooling stage, the temperature of the solidified droplet decreases due to the continuous cooling of the cold plate [83,94,95]. At this stage the ensuing droplet forms a peculiar shape with a pointed tip as shown in the last sequence of Figure 4.

**Figure 4.** Sequences of water droplet impact on a cold hydrophobic surface, in a way (**a**–**m**): the starting point of the process to the end of the phenomenon. Droplet diameter = 2.82 mm, impact velocity = 0.77 m/s, and surface temperature = −10 ◦C (reproduced with permission from [68]).

The deformation of an impacting water droplet throughout the freezing process on a hydrophobic horizontal surface is shown in Figure 4. Experimental studies have shown that sub-cooled temperature has no significant effect on the spreading of the droplet after impact [80,96]. However, it shows significant suppression of the receding speed and height when the period of receding is sufficiently long [96,97].

#### *2.2. E*ff*ect of Surface Wettability*

Surface wettability is the ability of surface to be wetted by a liquid and is mainly determined by surface roughness and chemistry. The Wenzel [98] and Cassie-Baxter [99] states describe the droplet wetting regime on the surface (see Figure 5). In the Wenzel state, water penetrates into the surface structures and conforms to the surface while in the Cassie-Baxter state, the water droplet remains above the surface structures, maintaining an almost spherical shape [4].

**Figure 5.** Wetting regimes of droplet on a solid surface: (**a**) Cassie-Baxter state and (**b**) Wenzel state (reproduced with permission from [4]).

Static contact angle, θ, which is defined as the angle between the droplet and the surface at the contact line (see Figure 6) is a characteristic used to define the degree of surface wettability. The advancing contact angle, θ*A*, is the contact angle for a droplet with an advancing contact line (e.g., for a growing droplet), while the receding contact angle, θ*R*, is the contact angle for a receding contact line (e.g., a shrinking droplet). Contact angle hysteresis is defined as θ*<sup>A</sup>* − θ*<sup>R</sup>* and is an indication of the droplet mobility on the surface, in a way that the lower is the hysteresis the easier the droplet moves on the surface (higher mobility). In general, in contact with water, a hydrophilic surface displays a contact angle of less than 90◦, while a hydrophobic surface shows a contact angle of more than 90◦. On a superhydrophobic surface, contact angle is more than 150◦, and the contact angle hysteresis should be less than 10◦ [4,68].

**Figure 6.** Schematic of static, θ, receding, θ*R*, and advancing, θ*A*, contact angles [100].

The wettability of a surface considerably affects the droplet dynamics after impact. It is shown that the wetted area during the spreading regime decreases with lower surface wettability [68]. It is also found that as wettability decreased a longer time was required for a droplet to reach equilibrium. This is due to oscillations in the droplet as more kinetic energy remains after the spreading regime. As spreading is suppressed, less energy is dissipated [68]. Tang et al. [77] studied the droplet dynamics of various liquid on surfaces with variable roughness. It was shown that as surface roughness increases, the droplet has a slower spreading time and smaller maximum spreading diameter. It was also observed that increase of surface roughness, promotes droplet splashing. Moreover, it was demonstrated that surface roughness has a prominent impact on the promotion of splashing in liquids with a smaller Ohnesorge number.

Surface wettability has a significant impact on the heat transfer rate in droplet dynamics. Pan et al. [68] studied droplet impacts on cold surfaces and observed how wettability affects the freezing process. As surface wettability decreases the total icing time of the droplet increases. This is because of a lower heat transfer rate due to a smaller heat transfer area. The lower surface contact area is due to the high interfacial tension in low wettability surfaces, which prevents droplet spreading. As expected, on hydrophilic surfaces the total freezing time is shorter comparing to the hydrophobic surfaces. On hydrophilic surfaces spreading is promoted due to more heat transfer area. Overall, more contact surface area leads to a higher heat transfer rate and in turn a shorter freezing time [68].

Surface wettability also has an impact on the boiling process during droplet dynamics. In general, as wettability increases the rate of phase-change heat transfer enhances [92,101]. As surface becomes more hydrophilic, smaller and faster bubble growth occurs [102,103]. Kim et al. [92] studied the effects of surface wettability on droplet rebounding on hot surfaces above the Leidenfrost temperature. In their experiments four surfaces were prepared, a smooth hydrophilic, a smooth hydrophobic, a hydrophilic, and a hydrophobic surface both with nanoscale structures. For both smooth and nanoscale hydrophobic surfaces, rebounding of the droplet was achieved, but at varying surface temperatures. For both hydrophilic surfaces, it was found that a higher surface temperature was needed to induce rebounding as more energy was required to disperse the droplet from the surface. On hydrophilic surfaces, there is more contact area and more surface tension, thus more work is needed to break the droplet adherence to the surface. Finally, on the nanoscale hydrophilic surface, splashing was induced caused by the capillary effects and cavities of the surface [92].

#### *2.3. E*ff*ect of Impact Velocity*

One of the most important characteristics affecting the dynamics of droplet impact is the droplet impact velocity. Specifically, impact velocity has a profound effect on the droplet spreading regime. Clearly, when the impact velocity increases, more kinetic energy is given to the droplet. The momentum in the vertical direction is then transferred in the radial direction upon impact [3,67,68]. As higher impact velocities a higher degree of spreading is promoted, leading to more surface area for heat transfer [62,66–68,93,104]. Rajesh et al. [67] showed that on superheated concave and convex surfaces, an increase in impact velocity translated to an increase in both maximum droplet spreading and contact time. This increase in maximum spreading diameter is due to a rise in impact kinetic energy allowing the droplet to spread more. Similarly, the increase in maximum droplet spreading, due to higher impact velocity, affects the freezing time of the droplet on cold surfaces. With higher impact velocities, a trend of shorter icing time and faster freezing is observed (faster spreading and higher heat transfer dispersal occurs upon impact [68]). Furthermore, as mentioned above, if the droplet has enough energy, it rebounds after retraction. Increase of the impact velocity has shown to promote rebounding and further increase of the impact velocity can induce splashing [3,59]. Chen et al. [93] studied the effect of the initial droplet impact momentum on the Leidenfrost point and detected that at high impact velocity, a higher surface temperature is needed to induce the Leidenfrost effect. The increased droplet impact momentum leads to higher kinetic energy in the spreading droplet, this causes a thinner liquid disk, hence nucleate bubbles are easier to burst out from the disk, which breaks the Leidenfrost state [93].

#### *2.4. E*ff*ects of Droplet Shape and Volume*

The shape of the droplet can impact its dynamics and the heat transfer rate. By increasing the ellipticity and asymmetry of the droplet, the spreading and retraction regimes are significantly changed. Generally, asymmetrical droplets require a higher impact velocity for rebounding [105–108]. Yun [109] studied the effects of asymmetry and ellipticity on droplet dynamics, specifically on its rebound regime. Asymmetric droplets, in the shape of an egg, showed uneven spreading and retraction in both the *x* and *y* axis (see Figure 7). To investigate the effect of asymmetry and ellipticity on the rebound regime, the height of the droplet at the center of mass was measured. It was found that with increasing ellipticity the height of the droplet decreased [109]. Additionally, with decreasing asymmetry the droplet showed a reduced rebound height. This suggests that suppression of droplet rebound can be induced by variances in its ellipticity and asymmetry. In other words, the degree of droplet deposition can be promoted with variances in the droplet's shape. The suppression of droplet bouncing is due to a break in the horizontal momentum of the droplet at low asymmetry and high ellipticity [109].

**Figure 7.** Simulation and experimental snapshots of the impact dynamics of (**a**) spherical and (**b**) asymmetric droplets on a hydrophobic surface (reproduced with permission from [109]); the equivalent droplet diameter is 2 mm and *We* = 17. SV and BV stand for the side and bottom-views in the experiment, respectively.

The size of the impacting droplet can also affect multiple characteristics of droplet dynamics. Pan et al. [68] studied how the initial droplet diameter affected its maximum spreading diameter and heat transfer rate. By increasing the initial droplet size, water droplet spreads more and showcases more oscillation during the spreading regime owing to an increase in initial kinetic energy. Additionally, upon impact on a cold surface, larger droplets showed increase in freezing time. Despite an increase in surface area during the spreading regime, larger droplets needed a longer time to freeze due to an increase in mass [68]. This experiment demonstrated the predominance of droplet volume over the effects of droplet dynamics on freezing time due to longer time spent during the heat transfer process

compared to the droplet dynamic process. It is worth mentioning that an increase in droplet volume can induce splashing due to an increase in the Weber number.

#### *2.5. E*ff*ect of Relative Humidity*

Relative humidity has considerable influence on droplet dynamics. Bobinski et al. [110] qualitatively observed the effects of relative humidity and surface temperature on droplet impact in icing conditions. In their studies, relative humidity was kept constant and equal to 75%, while the surface temperature was reduced to −10 ◦C. It was found that droplet freezing is triggered during the spreading phase and it freezes when the maximum spreading is reached, even on the superhydrophobic surface. By decreasing the relative humidity to 9%, it was shown that the impact process on the superhydrophobic surface is similar to the case of room temperature tests and the droplet does not freeze. Experimental study was also performed by Jadidi et al. [111] to understand the effects of relative humidity, surface temperature, and wettability, and the frosting mechanism on water droplet impact dynamics. In their study superhydrophobic, aluminum, and glass surfaces with three different surface temperatures (i.e., 20, 2, and −2 ◦C) were used. Furthermore, three different relative humidities (i.e., 10%, 20%, and 30%) were applied while the droplet Weber and Reynolds numbers, and the air temperature were fixed. It was revealed that the ratio of the surface temperature to the dew point temperature, which depends on relative humidity and air temperature, has a significant impact on droplet spreading, recoil, and contact angle. In general, when the mentioned ratio is less than one and decreases (it can be done by increasing the relative humidity or decreasing the surface temperature), condensation and frost formation become important, droplet spreading diameter increases significantly, and the equilibrium contact angle decreases.

#### *2.6. E*ff*ect of Dimensionless Numbers*

The effect of dimensionless numbers on the spreading ratio, *S*, during droplet deposition on a dry surface is shown in Figure 8 [60]. To obtain these results, one dimensionless number was individually changed while other dimensionless numbers as well as contact angles were kept constant. As shown in Figure 8a, the maximal spreading ratio increases with Reynolds number. On the other hand, the minimal spreading ratio at the end of the receding phase slightly decreases as the Reynolds number increases. Furthermore, the Reynolds number does not significantly affect the duration of spreading and receding phases in terms of dimensionless time. In contrast, as can be seen in Figure 8b, the Weber number does not affect the spreading phase. However, by increasing the Weber number, receding is slowed down. Increasing the Bond number, results in a slight increase of the spreading ratio in all impact stages (see Figure 8c). In short, the competition of inertial and viscous forces has significant influence on the spreading phase, while the surface tension force mainly affects the receding phase [60].

Herbert et al. [60] also showed that by increasing the Reynolds number, the total heat transfer from the substrate to the droplet significantly decreases. It was explained that although the heat transfer area is larger during the spreading and the beginning of the receding stages for large values of Reynolds, the heat flow is noticeably lower in these phases. In addition, it was shown that increasing the Weber number causes the dimensionless heat flow during the spreading and the receding phases to decrease and enhance, respectively. Moreover, increasing the Prandtl number results in significant reduction of heat transfer due to a suppressed convective heat transfer at the solid–liquid interface. It was also revealed that a higher value of Bond number causes the heat transfer rate to increase throughout the entire process since the wetted area is enlarged, as shown in Figure 8c.

**Figure 8.** Effect of dimensionless numbers on the droplet spreading ratio; (**a**) effect of the Reynolds number (*We* = 15, *Bo* = 0.5, *Pr* = 9.54), (**b**) effect of the Weber number (*Re* = 1000, *Bo* = 0.5, *Pr* = 9.54), and (**c**) effect of the Bond number (*We* = 15, *Re* = 1000, *Pr* = 9.54; reproduced with permission from [60]).

#### **3. Droplet Impact on Moving Surfaces**

Droplet impact onto moving surfaces shows vast differentiations in their dynamics. The degree of droplet deformation and induction of splashing are heavily dependent on droplet impact velocity and surface tangential velocity. Povarov et al. [112] observed the presence of the air layer causing the droplet to lift off at high tangential velocities. At low surface velocity, the droplet is deposited

onto the surface with no lift off. The bottom layer of the droplet adheres to the surface and spreads in a tear-like shape. At higher velocities, there is partial lift off of the droplet. The bottom layer of the droplet adheres to the surface and moves with it while the remainder of the droplet stays above the impact point. As the bottom layer is dragged away, sufficiently high tangential velocities cause the formation of an air layer. This triggers the rear edge of the droplet to splash and partially rebound off the surface. Total rebounding occurs at even higher tangential velocities as the droplet becomes significantly deformed [112]. Figure 9 illustrates droplet deformation on a moving surface. Mundo et al. [113] studied the splashing threshold of a droplet during impact on a moving substrate. Based on their observations, for high impact velocities, as the droplet impacts on the surface, a liquid film forms and spreads. The liquid film then deforms with the moving surface forming corona around the droplet. As the droplet continues to impact the surface, the fluid is fed into the corona film allowing it to grow. Whereas when less fluid is fed into the corona film, it becomes thinner and eventually disintegrates into secondary droplets. At lower impact velocities, the droplet impacts and spreads on the surface since there is not enough normal momentum to allow the corona formation and the kinetic energy is lost during the deformation process [86,113]. Mundo et al. also investigated the effects of surface roughness on droplet impacts on moving surfaces. They figured out that at increased surface roughness the droplet shows irregular deformation. Here the formation of the corona film and sequential instabilities were not observed due to vigorous splashing. Increased surface roughness also showed an increase in the formation of secondary droplets. On the other hand, decreased surface roughness promotes deposition of the droplet and splashing can be suppressed [113].

Zen et al. [86] studied how impacting velocity affected the splashing on moving surfaces. They described two splash phenomena, splash-around and side-splash in the rear edge. The threshold between splash-around and side-splash is dependent on the impact and surface velocities. In general, increasing the impact velocity induces side-splash and further increases in velocity lead to splash-around. It was also observed that increase of the impact velocity promotes the formation of smaller secondary droplets. Increasing the tangential velocity lead to promotion of side-splashing in the rear edge and lower tangential velocities promotes splash-around in all edges. Yao and Cai [114] demonstrated the

predominance of tangential velocity on splashing threshold. They observed that as the surface velocity increased, the critical impact velocity to induce splashing significantly decreased.

#### **4. Conclusions**

The importance of droplet impact dynamics in a vast number of industrial applications such as spray coating and painting, spray cooling, ink-jet printing, combustion engines, etc., drives numerous studies into this area. This paper summarized the intensive training course aimed at discussing some of the main governing factors on droplet impact dynamics and heat transfer on both static and moving solid surfaces. In details on solid static surfaces, the effects of surface temperature, wettability, impact velocity, dimensionless numbers, and droplet shape and volume were discussed. Surface temperature played a major role on droplet dynamics, particularly at super-heated or super-cooled temperatures. At ample temperatures, the Leidenfrost effect or the freezing process can be induced, wildly changing the droplet dynamics. Surface wettability also greatly affected the droplet dynamics as spreading, retraction, and rebounding can be promoted or suppressed depending on impact conditions. Impact velocity of the droplet heavily affects spreading and splashing regimes owing to transfer of kinetic energy after impact.

The effects of tangential and impact velocity on the impact dynamics on moving surfaces were also discussed in this article. Droplet impact dynamics are heavily dependent on the impact and surface velocity. Splashing is easily promoted on moving surfaces due to deformations in droplet caused by the moving surface. In general, having training on droplet impact subject for engineering students is essential due to the real-world applications of this topic.

**Author Contributions:** Conceptualization, supervision, writing-review and editing, S.M.; writing-original draft preparation, C.L. and M.J. All authors have read and agreed to the published version of the manuscript.

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

**Acknowledgments:** The authors gratefully acknowledge the support of the "Honors Program" at California State University Long Beach.

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

#### **References**


© 2020 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/).

## *Article* **A CFD Tutorial in Julia: Introduction to Laminar Boundary-Layer Theory**

**Furkan Oz \* and Kursat Kara**

School of Mechanical and Aerospace Engineering, Oklahoma State University, Stillwater, OK 74078, USA; kursat.kara@okstate.edu

**\*** Correspondence: foz@okstate.edu

**Abstract:** Numerical simulations of laminar boundary-layer equations are used to investigate the origins of skin-friction drag, flow separation, and aerodynamic heating concepts in advanced undergraduate- and graduate-level fluid dynamics/aerodynamics courses. A boundary-layer is a thin layer of fluid near a solid surface, and viscous effects dominate it. Students must understand the modeling of flow physics and implement numerical methods to conduct successful simulations. Writing computer codes to solve equations numerically is a critical part of the simulation process. Julia is a new programming language that is designed to combine performance and productivity. It is dynamic and fast. However, it is crucial to understand the capabilities of a new programming language before attempting to use it in a new project. In this paper, fundamental flow problems such as Blasius, Hiemenz, Homann, and Falkner-Skan flow equations are derived from scratch and numerically solved using the Julia language. We used the finite difference scheme to discretize the governing equations, employed the Thomas algorithm to solve the resulting linear system, and compared the results with the published data. In addition, we released the Julia codes in GitHub to shorten the learning curve for new users and discussed the advantages of Julia over other programming languages. We found that the Julia language has significant advantages in productivity over other coding languages. Interested readers may access the Julia codes on our GitHub page.

**Keywords:** CFD; Julia; Blasius; Hiemenz; Homann; Falkner–Skan; boundary-layer

#### **1. Introduction**

Computational fluid dynamics (CFD) simulation is one of the vital steps of the design of a product that includes fluid motion. Since fluid dynamics are extremely complex and the motion equations have nonlinear terms, usage of numerical approaches is inevitable to simulate or predict fluid motion. Predictions can also be done with experiments. However, it is usually costlier than a regular CFD simulation. On the other hand, fundamental knowledge about the flow, that is planned to simulate, is necessary because of the required numerical approach decision. Learning the canonical flows well is extremely important in this point because most of the complex flow consists of a combination of a couple of canonical flows. An airfoil CFD simulation can consist of boundary-layer flow over the smooth part of the airfoil, mixing layer flow where the tail ends, and blunt body flow in the wake region. One airfoil simulation consists of three different canonical flows. A full understanding of the canonical flows is extremely important to model an accurate airfoil simulation.

Prandtl [1] stated that some of the terms in the Navier–Stokes equations can be neglected for the boundary-layer flows. As a result of this assumption, well-known boundary-layer equations arose. These approaches are still valid after a century and the resultant system of equations inspired lots of researchers in their studies. One of them was Blasius who is the Ph.D. student of Prandtl. Blasius [2] worked on the same problem as Prandtl did. However, he aimed to overcome the enigma of turbulence by considering the phenomenon of boundary-layer flow explained by Hager [3]. He further

**Citation:** Oz, F.; Kara, K. A CFD Tutorial in Julia: Introduction to Laminar Boundary-Layer Theory. *Fluids* **2021**, *6*, 207. https://doi.org/ 10.3390/fluids6060207

Academic Editor: Ashwin Vaidya

Received: 13 May 2021 Accepted: 28 May 2021 Published: 3 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

simplified the boundary-layer equations for a flat plate. He assumed that the flow is parallel. In other words, the velocity component in the parallel direction is not zero and the velocity in the transverse direction is zero. Moreover, he represented the resultant system of equations with a third-order ordinary differential equation, which is known as the Blasius similarity solution.

The Falkner–Skan similarity solution is another laminar similarity solution. Cebeci [4] explains the Falkner–Skan equations named after V. M. Falkner and Sylvia W. Skan. Falkner and Skan generalized the Blasius similarity solution for non-parallel flows, such as wedge flows and corner flows. The resultant equation can be used to predict the boundary-layer thickness for wedge flow. It has to be noted that these assumptions are available in the laminar region. Once the turbulence occurs, both of these similarity solutions will be inaccurate. In the incompressible region, the Falkner–Skan equation can be used without additional equation; however, after the compressibility limit, the temperature effect must be introduced to the system as an additional equation. If the Falkner–Skan equation is solved with the energy equation, the boundary-layer profile can be obtained in the compressible region.

One of the other Ph.D. students of Prandtl is Karl Hiemenz who worked on Hiemenz flow which is a type of stagnation point flow. Hiemenz [5] formulated and calculated the stagnation point problem as explained in the Schlichting [6]. The problem was a special case of the Falkner–Skan similarity solution. Howarth [7] also worked on the same problem and concluded similar results. Another Ph.D. student of Prandtl, Fritz Homann [8] worked on the same problem for axisymmetric bodies as Schlichting [6] explains. Homann's similarity formulation was for a sphere while Hiemenz's similarity formulation was for a cylinder. Both similarity solutions are being widely used for stagnation flows.

The aforementioned papers are the origins of the boundary-layer theory. As it is mentioned, the boundary-layer theory is the origin of many problems, such as the laminar to turbulent boundary-layer transition and flow separation. The main concern of these researches is to increase the performance of the vehicle because the transition increases the heat transfer and the vehicle requires a better thermal protection system at high speeds due to increased heat transfer rate. On the other hand, flow separation may lead to a lift force loss on the wing as a result, the performance of aircraft decreases. Since the focus of the present paper is the fundamentals of laminar boundary-layer theory, the details will not be provided about advanced researches. However, readers who are interested in details may check [9–12] for subsonic boundary-layer transition, Ref. [13–21] for supersonic/hypersonic boundary-layer transition, and Ref. [22–24] for flow separation. The other researches where boundary-layer flow is involved are [25–30].

Understanding the aforementioned fundamental flows are crucial for a senior undergraduate student or a graduate student in order to simulate more complex flows. However, modeling these equations in a computer environment requires more than knowledge about canonical flows instead it requires knowledge about programming languages as well. There are several learning modules and papers [31–35] for computational fluid dynamics simulation coding; however, there is not enough publication for Julia language [36]. It is a relatively new coding language among the other coding languages such as Fortran, C/C++, Python, and MATLAB but it is getting popular fast because it is trying to fill the gap between the language of the state-of-art CFD codes, Fortran and C/C++, and straightforward/user-friendly languages, Python and MATLAB. Most of the state-of-art CFD codes are written in Fortran and C/C++ because it is so fast and random access memory usage (RAM) can be reduced drastically. On the other hand, MATLAB and Python have user-friendly and easy syntax which makes them favorites of students, with a price, which is the speed of computation. Julia is a language that tries to fill this gap between two different coding styles. It has a user-friendly environment, while also being fast. It provides a working space that helps users to write clear, high-level, generic code that resembles mathematical formulas. Julia's ability to combine high-performance with productivity makes it a great choice for researchers working in different areas.

In this paper, fundamental flow problems such as Hiemenz flow, Homann flow, Blasius flow, and Falkner–Skan flow will be derived from scratch and modeled in the Julia environment. The finite-difference discretization in space with Thomas algorithm as linear system solver is used. Outputs of the code are provided and they are compared with the literature. Additionally, the advantages of the Julia language over other languages will be discussed. The computer codes and the implementation instructions will help students to understand the fundamental flows which will provide insight for student's course work and researches. Using these examples, they can solve more complex flow types and develop their own codes. We make all these codes available on GitHub and they are accessible to everyone. We provide installation instruction for Julia and the required packages in Appendix A. The GitHub link of the codes also can be found in Appendix A.

#### **2. Laminar Boundary-Layer Theory**

Understanding canonical flows is crucial to understanding more complex flows such as flow over a wing or an airfoil. Hosseini et al. [37] studied flow over a wing section with a direct numerical simulation (DNS) study. They created a great video about how flow is formed over the aircraft. The video shows the development of the boundary-layer, how the trailing edge looks like a mixing-layer flow, and how Karman vortex street type of wake structures occurs. If the pre-study work is examined, it can be seen that the mesh structure is built on the flow prediction. Using a finer mesh on the critical regions is a key point of the CFD and it requires predictions about the possible flow behavior. The boundary-layer is the origin of many engineering problems in aerodynamics, including wing stall, the skin friction drag on an object, and the heat transfer that occurs in high-speed flight. In this present paper, boundary-layer theory will be examined under two subsections, which are laminar boundary-layer problems and stagnation point problems. In this present paper, fundamental fluid dynamics problems related to boundary-layer theory will be derived and implemented in a relatively new programming language, Julia. The contribution is employing the Julia language. The governing equations and solution methodologies are already published in the literature. The interested reader should refer to the additional references [6,38] for detailed derivations. The manuscript may enable students to adopt the programming language easily.

#### *2.1. Laminar Boundary-Layer Flow Problems*

Velocity distribution over a flat plate can be represented with a similarity solution. In this subsection, Blasius and Falkner–Skan similarity solutions will be derived from scratch and they will be solved numerically in the Julia environment. The Julia codes will be available to shorten the learning curve.

#### 2.1.1. Blasius Flow Problem

Schematic description of flow over a flat plate, in other words, Blasius boundarylayer flow can be illustrated as in Figure 1. The governing equations for Blasius flow are boundary-layer equations. These equations can be obtained by non-dimensionalizing the two-dimensional Navier–Stokes equations that are:

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \tag{1}$$

$$
\mu \frac{\partial \mu}{\partial \mathbf{x}} + v \frac{\partial \mu}{\partial y} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}} + \nu \left( \frac{\partial^2 \mu}{\partial \mathbf{x}^2} + \frac{\partial^2 \mu}{\partial y^2} \right) \tag{2}
$$

$$
\mu \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = -\frac{1}{\rho} \frac{\partial p}{\partial y} + \nu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right), \tag{3}
$$

where *u* and *v* are the velocities in the *x*-direction and *y*-direction, respectively. *p* is the pressure, *ν* is the kinematic viscosity, *ρ* is the density. If Equations (1)–(3) are nondimensionalized with:

$$u^\* = \frac{u}{lL\_{\infty}}, \quad v^\* = \frac{v}{lL\_{\infty}}, \quad p^\* = \frac{p}{\rho Ul\_{\infty}^2}.\tag{4}$$

$$\mathbf{x}^\* = \frac{\mathbf{x}}{L}, \quad \mathbf{y}^\* = \frac{\mathbf{y}}{L}, \tag{5}$$

where *L* is the plate length and *U*∞ is the free-stream velocity. The two-dimensional, incompressible non-dimensional Navier–Stokes equations can be shown as:

$$
\frac{\partial u^\*}{\partial x^\*} + \frac{\partial v^\*}{\partial y^\*} = 0 \tag{6}
$$

$$u^\* \frac{\partial u^\*}{\partial x^\*} + v^\* \frac{\partial u^\*}{\partial y^\*} = -\frac{\partial p^\*}{\partial x^\*} + \frac{1}{Re\_{\infty}} \left( \frac{\partial^2 u^\*}{\partial x^{\*\frac{\nu}{2}}} + \frac{\partial^2 u^\*}{\partial y^{\*\frac{\nu}{2}}} \right) \tag{7}$$

$$u^\* \frac{\partial v^\*}{\partial x^\*} + v^\* \frac{\partial v^\*}{\partial y^\*} = -\frac{\partial p^\*}{\partial y^\*} + \frac{1}{Re\_{\infty}} \left( \frac{\partial^2 v^\*}{\partial x^{\*\prime}} + \frac{\partial^2 v^\*}{\partial y^{\*\prime}} \right). \tag{8}$$

Star superscript ([ ]∗) corresponds to non-dimensional variable. In a boundary-layer flow, some variables are smaller than others. For example, *u*∗ = *O*(1) and *x*∗ = *O*(1) so *∂u*∗ *<sup>∂</sup>x*<sup>∗</sup> = *O*(1). In the continuity equation, two terms must be in the same order to obtain 0 so *<sup>∂</sup>v*<sup>∗</sup> *<sup>∂</sup>y*<sup>∗</sup> = *O*(1). It can be seen that from Figure 1, *y* = *O*(*δ*∗) where *δ*<sup>∗</sup> is the nondimensionalized boundary-layer thickness. From here, it can be concluded that *v* = *O*(*δ*∗). It has to be noted that *δ*∗ << 1. If the same approach is applied to momentum equations, the final system of equations will be:

$$
\frac{
\partial 
\mu^\*
}{
\partial x^\*
} + \frac{
\partial v^\*
}{
\partial y^\*
} = 0
\tag{9}
$$

$$
\mu^\* \frac{\partial \mu^\*}{\partial x^\*} + \upsilon^\* \frac{\partial \mu^\*}{\partial y^\*} = -\frac{\partial p^\*}{\partial x^\*} + \frac{1}{\text{Re}\varepsilon\_{\infty}} \frac{\partial^2 \mu^\*}{\partial y^{\*2}} \tag{10}
$$

$$\frac{\partial p^\*}{\partial y^\*} = 0.\tag{11}$$

**Figure 1.** Schematic description of the flow over a flat plate.

For an inviscid flow, where the Reynolds number is so high, the viscous term can be neglected. The pressure is constant as obtained in the Equation (11) and *v* = 0. Equation (10) will be:

$$-\frac{dP\_\varepsilon^\*}{d\mathbf{x}^\*} = \mathcal{U}\_\varepsilon^\* \frac{d\mathcal{U}\_\varepsilon^\*}{d\mathbf{x}^\*}.\tag{12}$$

where *Ue* and *Pe* are the inviscid velocity and pressure. Finally, the two-dimensional boundary-layer equations in the dimensionless form can be written as:

$$
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \tag{13}
$$

$$
\mu \frac{\partial \mu}{\partial x} + \upsilon \frac{\partial \mu}{\partial y} = \mathcal{U}\_{\mathfrak{e}} \frac{d \mathcal{U}\_{\mathfrak{e}}}{\partial x} + \upsilon \frac{\partial^2 \mu}{\partial y^2}. \tag{14}
$$

It has to be noted that star superscript is suppressed to avoid confusion. In order to solve this problem, similarity transformation can be used. The similarity parameter (*η*) can be chosen as:

$$
\eta = \frac{y}{\sqrt{\frac{\nu \chi}{\Pi\_{\infty}}}}.\tag{15}
$$

After that, velocity components can be shown as:

$$
\mu(\mathbf{x}, \mathbf{y}) = \mathcal{U}\_{\infty} \mathcal{F}(\boldsymbol{\eta}) \tag{16}
$$

$$w(x, y) = \frac{G(\eta)}{\sqrt{\frac{x}{\Pi\_{\infty}\nu}}}.\tag{17}$$

If Equations (13) and (14) are rearranged with the given velocities, the boundary-layer equations will be:

$$G' = \frac{\eta}{2} F'\tag{18}$$

$$F'' = -\frac{1}{2}\eta F'' + GF'.\tag{19}$$

In order to make the next procedure easier, it can be assumed that *F* = *f* - . After that, Equation (18) can be integrated and put into Equation (19). The final ordinary differential equation can be obtained as:

$$2f'''' + ff'' = 0.\tag{20}$$

The boundary conditions of the system will be:

$$
\mu(\mathbf{x},0) = \mathcal{U}\_{\infty}F(\eta=0) = 0 \longrightarrow F(0) = f'(0) = 0 \tag{21}
$$

$$v(\mathbf{x},0) = \frac{G(\eta=0)}{\sqrt{\frac{\mathbf{x}}{\Pi\_{\text{av}}\mathbf{v}}}} = 0 \longrightarrow G(0) = f(0) = 0 \tag{22}$$

$$\mu(\mathbf{x}, \mathbf{y} \to \infty) = \mathcal{U}\_{\infty}^{\cdot} F(\eta \to \infty) = 1 \longrightarrow F(\infty) = f'(\infty) = f(\infty) = 1. \tag{23}$$

The given assumptions reduced the system of second-order partial differential equations into a third-order ordinary differential equation. The computational approach will be examined in the Julia framework.

#### 2.1.2. Numerical Solution of Blasius Flow Problem

In this computational section, the solution of the third-order ordinary differential equation will be transformed into two equations and they will be solved with Thomas algorithm [39] which is one of the best methods for the tridiagonal matrices. It has to be noted that the same equations can be solved with any other methods such as Runge-Kutta, Runge-Kutta-Fehlberg, compact finite difference, high-order finite-difference; however, in this present paper, authors used the finite difference method for spatial discretization with the Thomas algorithm for the linear system solution. If it is assumed that *f* - = *h* and *f* = *p* in Equation (20), the system of equations will be:

$$2\mathbf{h}'''' + p\mathbf{h}' = 0\tag{24}$$

$$p'-h=0,\tag{25}$$

with the following boundary conditions:

$$p(0) = 0\tag{26}$$

$$h(0) = 0\tag{27}$$

$$h(\infty) = 1.\tag{28}$$

If the second-order central finite difference method is applied to *h* variable:

$$h\_{n}^{\prime\prime} = \frac{h\_{n+1} - 2h\_n + h\_{n-1}}{(\Delta \eta)^2} + O((\Delta \eta)^2) \tag{29}$$

$$h\_n' = \frac{h\_{n+1} - h\_{n-1}}{2\Delta\eta} + O((\Delta\eta)^2). \tag{30}$$

In this paper, details of the finite difference and derivation of it from Taylor's series will not be covered. If the reader is curious about the derivation of them, the book of Moin [40] can be checked. If the finite difference approach for *h* is substituted into the Equation (24), the final equation will be:

$$A\_{\rm nl}h\_{\rm n+1} + B\_{\rm nl}h\_{\rm n} + C\_{\rm nl}h\_{\rm n-1} = 0,\tag{31}$$

where the *An*, *Bn*, and *Cn* are:

$$A\_n = \frac{2}{\Delta \eta^2} + \frac{p\_n}{2\Delta \eta} \tag{32}$$

$$B\_{ll} = \frac{-4}{\Delta \eta^2} \tag{33}$$

$$C\_n = \frac{2}{\Delta \eta^2} - \frac{p\_n}{2\Delta \eta},\tag{34}$$

with the boundary conditions for the system the tridiagonal system (*Ax* = *b*) can be shown as:

$$
\underbrace{\begin{bmatrix}1&0&0&0&\cdots&0\\A\_{2}&B\_{2}&C\_{2}&0&\cdots&0\\\vdots&\ddots&\ddots&\ddots&\vdots\\0&0&A\_{N-2}&B\_{N-2}&C\_{N-2}&0\\0&0&0&A\_{N-1}&B\_{N-1}&C\_{N-1}\\0&0&0&0&0&1\end{bmatrix}}\_{A}\underbrace{\begin{bmatrix}h\_{1}\\h\_{2}\\\vdots\\h\_{N-2}\\h\_{N-1}\\h\_{N}\end{bmatrix}}\_{x}=\underbrace{\begin{bmatrix}0\\0\\0\\\vdots\\0\\0\\1\end{bmatrix}}\_{b}\tag{35}
$$

It has to be noted that the first and last rows of the matrix are known and they are coming from the boundary conditions. In order to solve this problem, the Thomas algorithm can be applied. If relation between *hn* and *hn*−<sup>1</sup> is:

$$h\_n = G\_n + H\_n h\_{n+1} \tag{36}$$

$$h\_{n-1} = G\_{n-1} + H\_{n-1} h\_n. \tag{37}$$

If Equation (37) is substituted into Equation (31), the final equation will be:

$$h\_{\rm nl} = \frac{-\mathbb{C}\_{\rm n}\mathbb{G}\_{\rm n-1}}{B\_{\rm n} + \mathbb{C}\_{\rm n}H\_{\rm n-1}} + \frac{-A\_{\rm n}}{B\_{\rm n} + \mathbb{C}\_{\rm n}H\_{\rm n-1}}h\_{\rm n+1}.\tag{38}$$

If Equations (36) and (38) are compared, *Gn* and *Hn* coefficients can be found as:

$$G\_{\rm nl} = \frac{-\mathbb{C}\_{\rm n} G\_{\rm n-1}}{B\_{\rm n} + \mathbb{C}\_{\rm n} H\_{\rm n-1}} \tag{39}$$

$$H\_{\rm II} = \frac{-A\_{\rm II}}{B\_{\rm II} + \mathbb{C}\_{n}H\_{n-1}}.\tag{40}$$

The problem here is that one has to know *G*<sup>1</sup> and *H*<sup>1</sup> in order to start the calculation. *H*<sup>1</sup> can be assumed as 0 and the *G*<sup>1</sup> can be calculated as 0 from the boundary conditions. Once *G*<sup>1</sup> and *H*<sup>1</sup> are calculated, *Gn* and *Hn* can be calculated from Equations (39) and (40). The given formulations can be implemented as it is shown in Listing 1.

The problem in Listing 1 is the usage of the *p* and *h* values which are not defined yet. In order to start to iteration, the initial assumption for *h* must be given. Most of the time, the linear assumption is the best assumption. Figure 2 shows the schematic of the profile transformation after some iterations. Once the initialization of *h* is done, *p* can be calculated from the integration of *p*-= *h*. The discrete integration of *p* can be written as:

$$p\_n = p\_{n-1} + \int\_{\eta\_{n-1}}^{\eta\_n} h d\eta. \tag{41}$$

```
Listing 1. Implementation of Thomas algorithm for Blasius profile. ✞ ☎
```

```
1 A = [ 2/Δη2 + p[i]/(2*Δη) for i=1:N]
2 B=[−4/Δη2 for i=1:N]
3 C = [ 2/Δη2 − p[i]/(2*Δη) for i=1:N]
4 D=[0 for i=1:N]
5
6 for i=2:N−1
7 G[i] = − ( C[i]*G[i−1] + D[i] )/(B[i] + C[i] * H[i−1])
8 H[i] = − A[i] /(B[i] + C[i] * H[i−1])
9 end
10
11 for i=N-1:-1:2
12 h[i] = G[i] + H[i] * h[i+1]
13 end
```
**Figure 2.** Schematic description of the change in the velocity profile from initial guess to final solution.

✝ ✆

If the integration is done numerically with the trapezoidal rule [41], the initial guess for *h* and *p* can be calculated in Julia as shown in Listing 2

After the initialization is done, Listing 1 can be run with *p* calculation of Listing 2, until the change in *<sup>h</sup>* is smaller than an arbitrary parameter which can be taken as 1 <sup>×</sup> <sup>10</sup><sup>−</sup>8. The final profile will be as shown in Figure 3. The profile is also compared with Schlichting [6] in order to validate the results. It has to be noted that, Schlichting used *η* = *<sup>y</sup>* 2*νx U*∞ as

✝ ✆

similarity coordinate so the figure is plotted according to this nondimensionalization.

**Listing 2.** Implementation of initial velocity guess and initial *<sup>p</sup>* for Blasius profile. ✞ ☎

```
1 h = [(i−1)/(N−1) for i=1:N]
2
3 for i = 2:N
4 p[i] = p[i−1] + (h[i] + h[i−1])*Δη/2
5 end
```
**Figure 3.** The velocity distribution of the Blasius similarity solution obtained by given Julia code and data digitized from Schlichting [6].

#### 2.1.3. Falkner–Skan Flow Problem

The Falkner–Skan similarity solution can be considered as a family of the similarity solutions since it is a general solution that includes Blasius flow, Hiemenz flow (see Section 2.2.1), and more. Falkner–Skan equation cannot represent Homann flow (see Section 2.2.3) because Homann flow is an axisymmetric flow. On the other hand, Falkner– Skan flow is a two-dimensional flow. The name 'similarity solution' arises from the solutions at two arbitrary stations which are related to one another by means of a scale factor. A flow that can be represented with boundary-layer Equations (13) and (14), and satisfies the following equation:

$$\frac{\mu(\mathbf{x\_1})\_\prime \frac{y}{\mathcal{J}(\mathbf{x\_1})}}{\mathcal{U}\_\mathbf{f}(\mathbf{x\_1})} = \frac{\mu(\mathbf{x\_2})\_\prime \frac{y}{\mathcal{J}(\mathbf{x\_2})}}{\mathcal{U}\_\mathbf{f}(\mathbf{x\_2})} \, , \tag{42}$$

can be considered as self-similar. If boundary-layer Equations (13) and (14) are considered and similarity a transformation is assumed as:

$$
\mu(\mathbf{x}, \mathbf{y}) = \mathcal{U}\_{\mathbf{c}}(\mathbf{x}) f'(\mathbf{f}), \ \mathbf{y} = \frac{\mathcal{Y}}{\mathcal{g}(\mathbf{x})}.\tag{43}
$$

If the continuity equation of boundary-layer Equation (13) is modified with the given transformation, the final equation will be:

$$\frac{\partial \upsilon}{\partial \xi} = -\mathsf{U}\_{\mathsf{e}}' f' \mathsf{g} + \mathsf{U}\_{\mathsf{e}} f'' \mathsf{g}' \xi. \tag{44}$$

If Equation (44) is integrated over *ξ* to find the velocity *v*, the *v* velocity will be:

$$
\sigma(\mathbf{x}, \mathfrak{F}) = -(\mathcal{U}\_{\mathfrak{e}} \mathfrak{g})' f + \mathcal{U}\_{\mathfrak{e}} \mathfrak{F} \mathfrak{g}' f' + H(\mathbf{x}).\tag{45}
$$

*H*(*x*) can be calculated from the boundary conditions on the wall where *ξ* is zero and correspondingly, velocities and *f* - (0) are zero. *H*(*x*) will be 0 if *f*(0) is chosen as 0. If calculated derivatives and velocities are substituted into Equation (14), after some simplification, the final equation will be:

$$f^{\prime\prime\prime} + \frac{\mathcal{g}}{\nu} (\mathcal{U}\_c \mathcal{g})^{\prime} f f^{\prime\prime} + \frac{\mathcal{g}^2 \mathcal{U}\_c^{\prime}}{\nu} (1 - f^{\prime}) = 0. \tag{46}$$

Since *f* must be the function of *ξ* only and must not be a function of *x*, the coefficients of second and third terms must be constant. The final Falkner–Skan equation for the family of self-similar solutions:

$$f^{\prime\prime\prime} + a f f^{\prime\prime} + \beta (1 - f^{\prime 2}) = 0,\tag{47}$$

where *α* = *<sup>g</sup> <sup>ν</sup>* (*Ueg*) and *<sup>β</sup>* <sup>=</sup> *<sup>g</sup>*2*U*- *e <sup>ν</sup>* . The *α* and *β* can be further solved for the velocity scale and length scale. If the following consideration is applied:

$$\mathfrak{n}^\* = \mathfrak{g}(\mathcal{U}\_\mathfrak{e}\mathfrak{g})'\tag{48}$$

$$
\mathcal{S}^\* = \mathcal{g}^2 \mathcal{U}\_c'. \tag{49}
$$

If (*Ueg*2) is written in terms of *α*∗ and *β*∗, the obtained equation integrated with respect to *x*, the resultant equation will be:

$$\text{dL}\_{\text{f}}\text{g}^{2} = (\text{2a}^{\*} - \beta^{\*})\text{x} + c.\tag{50}$$

The constant of the integration represents a shift in the origin on x. Hence it doesn't affect the result and it also can be calculated from the stagnation point where *x* = 0 and *Ue* = 0 as a result, *c* = 0. If relation of *β*<sup>∗</sup> = *g*2*U*- *<sup>e</sup>* is divided by Equation (50) and integrated with respect to *x*, *Ue* can be calculated as:

$$
\mathcal{U}\_{\varepsilon} = \varepsilon\_1 \mathbf{x}^{\mathfrak{m}},
\tag{51}
$$

where *m* = *<sup>β</sup>*<sup>∗</sup> <sup>2</sup>*α*∗−*β*<sup>∗</sup> and *<sup>c</sup>*<sup>1</sup> is a positive or negative constant which depends on the sign of *Ue*. It can be concluded from these calculations that similar solutions exist when the inviscid velocity is proportional to *x* raised to some power. Next, Equation (50) can be used by taking *c* = 0 to calculate the *g* which is:

$$\mathbf{g} = \sqrt{\frac{2a^\*}{c\_1(1+m)}} \mathbf{x}^{1-m}.\tag{52}$$

Self-similar boundary-layers occur when the external velocity is the simple power law (*Ue* = *U*0(*x*/*L*)*m*), where the arbitrary constants *U*<sup>0</sup> and L have the same sign as *U* and x. The similarity variable for these kinds of flows can be written as:

$$\eta = \frac{y}{\delta} = \frac{y}{\sqrt{\pm \frac{\nu \chi}{\Pi\_c}}} = \frac{y}{\sqrt{\pm \frac{\nu L}{\Pi\_0} \left(\frac{x}{L}\right)^{1-m}}}.\tag{53}$$

When *Ue* and *x* have the same signs, the Falkner–Skan equation can be written as:

$$f^{\prime\prime\prime} + \frac{1}{2}(m+1)f f^{\prime\prime} + m(1 - f^{\prime2}) = 0.\tag{54}$$

The two arbitrary constants *α* and *β* have been reduced to one constant *m* by fixing the scale for the function *δ*(*x*). The boundary conditions of the equation are:

$$f(0) = 0\tag{55}$$

$$f'(0) = 0\tag{56}$$

$$f'(\eta \to \infty) = 1.\tag{57}$$

If the Falkner–Skan Equation (54) is carefully examined, it can be seen that when the constant *m* is 0, the equation will be Blasius flow. Moreover, if the constant *m* is 1, the equation will be Hiemenz flow (see Section 2.2.1). This important point is stated before and it is emphasized one more time after the derivation. Another great usage of the Falkner– Skan equation is to simulate the boundary-layer over a wedge with half-angle *<sup>θ</sup>* <sup>=</sup> *<sup>m</sup><sup>π</sup> <sup>m</sup>*+<sup>1</sup> when the *m* is between 0 and 1. If the *m* is in between 1 and 2, the Falkner–Skan equation

will solve a corner flow with *θ* > *<sup>π</sup>* <sup>2</sup> . The visual schematic of four different physical flows that can be calculated from the Falkner–Skan equation can be seen in Figure 4. One interesting point of the equation is that when *m* = −0.0904, the obtained profile will have zero-shear at the wall which corresponds to the verge of the separation point for all *x* stations.

**Figure 4.** The representation of flow types that can be calculated with the Falkner–Skan equation. (**a**) Blasius flow (*m* = 0), (**b**) wedge flow (0 < *m* < 1), (**c**) Hiemenz flow (*m* = 1), and (**d**) corner flow (1 < *m* < 2).

2.1.4. Numerical Solution of Falkner–Skan Flow Problem

In this section, the Falkner–Skan equation will be solved with the Thomas algorithm and central finite difference scheme. In order to do that, the third-order ordinary differential equation should be reduced to second-order and first-order differential equations as it is done in Blasius flow. If it is assumed that *f* - = *h* and *f* = *p* in Equation (54), the system of equations will be:

$$h'' + \frac{1}{2}(m+1)ph' + m(1-h^2) = 0\tag{58}$$

$$p'-h=0,\tag{59}$$

where the boundary conditions of the system are:

$$v(0) = 0\tag{60}$$

*h*(0) = 0 (61)

$$h(\infty) = 1\tag{62}$$

If finite difference scheme is applied to h variable, the final system of equations will be:

$$A\_n h\_{n+1} + B\_n h\_n + C\_n h\_{n-1} + D\_n = 0,\tag{63}$$

where the *An*, *Bn*, *Cn*, and *Dn* are:

$$A\_n = \frac{1}{\Delta \eta^2} + (m+1)\frac{p\_n}{4\Delta \eta} \tag{64}$$

$$B\_n = \frac{-2}{\Delta \eta^2} - mh\_n \tag{65}$$

$$C\_n = \frac{1}{\Delta \eta^2} - (m+1)\frac{p\_n}{4\Delta \eta} \tag{66}$$

$$D\_{\rm nl} = m.\tag{67}$$

The relation between *hn* and *hn*−<sup>1</sup> is taken as it is done in Blasius flow. When Equation (37) is substituted into Equation (63), the final *Gn* and *Hn* will be:

$$G\_n = \frac{-\mathbb{C}\_n G\_{n-1} + D\_n}{B\_n + \mathbb{C}\_n H\_{n-1}}\tag{68}$$

$$H\_{\rm II} = \frac{-A\_{\rm II}}{B\_{\rm II} + C\_{\rm II}H\_{\rm II-1}}.\tag{69}$$

The implementation of these variables in the Julia environment can be seen in Listing 3.

**Listing 3.** Implementation of the Thomas algorithm for the Falkner–Skan profile. ✞ ☎

```
1 m = 0.5
2 h = [(i−1)/(N−1) for i=1:N]
3
4 for i = 2:N
5 p[i] = p[i−1] + (h[i] + h[i−1])*Δη/2
6 end
7
8 while 1e-8<=errorProfile
9
10 A = [ 1/Δη2 + (m+1)*p[i]/(4*Δη) for i=1:N]
11 B=[−2/Δη2 − m*h[i] for i=1:N]
12 C = [ 1/Δη2 − (m+1)*p[i]/(4*Δη) for i=1:N]
13 D=[m for i=1:N]
14
15 for i=2:N-1
16 G[i] = − ( C[i]*G[i−1] + D[i] )/(B[i] + C[i] * H[i−1])
17 H[i] = − A[i] /(B[i] + C[i] * H[i−1])
18 end
19
20 hp = copy(h)
21
22 for i=N−1:−1:2
23 h[i] = G[i] + H[i] * h[i+1]
24 end
25
26 errorProfile = maximum( abs .(hp -h))
27
28 for i = 2:N
29 p[i] = p[i−1] + (h[i] + h[i−1])*Δη/2
30 end
31 end
```
The final profiles for the varying *m* values can be seen in the Figure 5. The profiles are also compared with Schlichting [6] in order to validate the results. It has to be noted that, Schlichting used *η* = *m*+1 <sup>2</sup> √ *y νx U*∞ as similarity coordinate so the figure is plotted according to this nondimensionalization.

**Figure 5.** The velocity distribution of the Falkner–Skan similarity solution obtained by Julia code and data digitized from Schlichting [6] for varying *m* values.

#### *2.2. Stagnation Point Flow Problems*

Boundary-layer velocity distribution over a wall that is perpendicular to the flow velocity vector can be represented with a similarity solution. In this subsection, Hiemenz and Homann similarity solutions will be derived from scratch and they will be solved numerically in the Julia environment. The Julia codes will be available to shorten the learning curve.

#### 2.2.1. Hiemenz Flow Problem

This problem is that of a fluid flow that is parallel to the *y*-axis in the far-field impinging on a wall that coincides with the *x*-axis. The flow which is perpendicular to a cylinder can be assumed as the Hiemenz flow around the stagnation point. The schematic description of the Hiemenz flow can be seen in Figure 6. In the Hiemenz flow, viscous forces away from the wall become so small in comparison with the inertia forces, particularly when the Reynolds number is large. In this case, inviscid irrotational flow assumption (*ξ* = ∇ × *U* = 0) can be done. The velocity can be represented with a scalar function, *φ*, which is the velocity potential. The velocities in the *x*-direction and *y*-direction can be written as:

**Figure 6.** Schematic description of the Hiemenz and Homann flow. **Left** side corresponds to general flow and the **Right** side is the extended vision of the dashed rectangle. The flow represents Hiemenz flow if the circle is the projection of a cylinder and Homann flow if the circle is the projection of a sphere.

$$
\mu = \frac{\partial \phi}{\partial x}, \ v = \frac{\partial \phi}{\partial y}. \tag{70}
$$

If the dimensional Navier–Stokes equations (Equations (1)–(3)) are considered for this flow as well, and the continuity equation becomes:

$$
\frac{
\partial^2 \phi
}{
\partial x^2
} + \frac{
\partial^2 \phi
}{
\partial y^2
} = 0.
\tag{71}
$$

On the wall, in the absence of viscosity, the flow can slip in the *x*-direction but in *y*-direction the velocity must be zero ( *∂φ <sup>∂</sup><sup>y</sup>* |*y*=<sup>0</sup> = 0) because of the no-penetration boundary condition. The inviscid flow solution for the potential function was found to be:

$$\phi\_{(x,y)} = \frac{a}{2}(x^2 - y^2),\tag{72}$$

where *a* is a constant that depends on the freestream flow and the body shape. This solution satisfies the governing equations for inviscid, irrotational flow, and the boundary conditions which can be shown as:

$$
\mu = \frac{\partial \phi}{\partial x} = a\mathbf{x}, \ v = \frac{\partial \phi}{\partial y} = -ay. \tag{73}
$$

Since potential the function is obtained, the stream-function can be calculated from the potential function. the velocity components, in terms of stream-function, *ψ*(*x*,*y*):

$$u = a\mathbf{x} = \frac{\partial \psi}{\partial y'}\quad v = -ay = -\frac{\partial \psi}{\partial \mathbf{x}}.\tag{74}$$

If Equation (74) is integrated for *x* and *y*, the stream function can be found as:

$$
\psi = axy + c.\tag{75}
$$

In order to find the pressure in the inviscid flow, Bernoulli equation [42] can be used. The pressure from the Bernoulli equation is:

$$p\_0 - p = \frac{a^2 \rho}{2} (x^2 + y^2),\tag{76}$$

where *p*<sup>0</sup> is the stagnation pressure and *ρ* is the density. So far, flow is assumed as inviscid flow; however viscous forces will modify the inviscid solution, in particular in the *y*direction. Hence, the viscous velocity components can be assumed as:

$$
\mu(x, y) = \text{xg}(y) \tag{77}
$$

$$v(x, y) = -f(y) \tag{78}$$

If Equations (77) and (78) are substituted in the Navier–Stokes equations (Equations (1)–(3)), the final system of equations will be:

$$\mathbf{g} = \mathbf{f}'\tag{79}$$

$$\text{xg}^2 - \text{xfg}' = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}} + \nu \text{xg}'' \tag{80}$$

$$ff' = -\frac{1}{\rho} \frac{\partial p}{\partial y} + \nu f''. \tag{81}$$

Once Equation (79) is substituted into Equations (80) and (81), two non-linear ordinary differential equations can be obtained as:

$$\|\mathbf{x}f'^2 - \mathbf{x}f\mathbf{y}\prime = -\frac{1}{\rho}\frac{\partial p}{\partial \mathbf{x}} + \nu \mathbf{x}f\prime\prime\tag{82}$$

$$ff' = -\frac{1}{\rho} \frac{\partial p}{\partial y} + \upsilon f'' \tag{83}$$

It has to be noted that these two equations are coupled, which means they have to be solved together. If the *x*-momentum and *y*-momentum are solved for pressure, it can be seen that the pressure will be in the form of:

$$p(x,y) = \frac{1}{2}x^2(constant) + H(y),\tag{84}$$

where *H*(*y*) is a function depends on only *y*. On the other hand, if the inviscid pressure compared with the viscous pressure as *y* −→ ∞, the pressure gradient in *x*-direction and *y*-direction can be found as:

$$\frac{\partial p}{\partial \mathbf{x}} = -a^2 \rho \mathbf{x} \tag{85}$$

$$\frac{\partial p}{\partial y} = -\frac{a^2 \rho}{2} F'(y). \tag{86}$$

If the pressure gradients are substituted into Equations (82) and (83), the resultant system of ordinary differential equations can be found as:

$$f'^2 - ff'' = a^2 + \nu f''''\tag{87}$$

$$ff' = \frac{a^2}{2}F' - \nu f''.\tag{88}$$

Momentum equations are decoupled in this system of equations so once, *f*(*y*) is calculated from Equation (87), *F*(*y*) in the Equation (88) can be solved. The boundary conditions of the system are:

$$
\mu(x,0) = 0 \to f'\_{(y=0)} = 0 \tag{89}
$$

$$
\sigma(x,0) = 0 \to f\_{(y=0)} = 0 \tag{90}
$$

$$
\mu(\mathbf{x}, \mathbf{y} \to \infty) = a\mathbf{x} \to f\_{\infty}' = a \tag{91}
$$

$$v(x, y \to \infty) = -ay \to f\_{\infty}^{\prime\prime} = 0\tag{92}$$

Additionally, *<sup>F</sup>*(*<sup>y</sup>* −→ <sup>∞</sup>) = *<sup>y</sup>*<sup>2</sup> can be obtained from *<sup>y</sup>*-momentum equation as *<sup>y</sup>* −→ <sup>∞</sup>. The final Equations (87) and (88) can be solved with the given boundary conditions; however, it is possible to get rid of the dependence on *a*. In order to do that, affine transformation [43] can be used. If it is assumed that *f*(*y*) = *Aϕ*(*η*) and *y* = *Bη* where *A* and *B* are constant to be determined, the *f* derivative functions can be calculated as:

$$f' = \frac{A}{B} \varphi'\tag{93}$$

$$f^{\prime\prime} = \frac{A}{B^2} \varphi^{\prime\prime} \tag{94}$$

$$f^{\prime\prime\prime} = \frac{A}{B^3} \varphi^{\prime\prime\prime} \tag{95}$$

If Equations (93)–(95) substituted into the Equations (87) and (88) and assume that *ν AB* <sup>=</sup> 1 and *<sup>a</sup>*2*B*<sup>2</sup> *<sup>A</sup>*<sup>2</sup> = 1, the *<sup>A</sup>* and *<sup>B</sup>* coefficients can be calculated as:

$$A = \sqrt{av} \tag{96}$$

$$B = \sqrt{\frac{\nu}{a}}\tag{97}$$

The given transformation allows to reduce two ordinary differential equations into one third-order ordinary differential equation as:

$$
\varphi^{\prime\prime\prime} + \varphi \varphi^{\prime\prime} - \varphi^{\prime 2} + 1 = 0 \tag{98}
$$

where the boundary conditions are:

$$f'(0) = 0 \to \varphi'(0) = 0\tag{99}$$

$$f(0) = 0 \to \varphi(0) = 0\tag{100}$$

$$f'(\infty) = a \to \varphi'(\infty) = 1\tag{101}$$

and the velocities:

$$
\mu = a \mathbf{x} \boldsymbol{\varphi}'(\boldsymbol{\eta}) \tag{102}
$$

$$w = -\sqrt{a\nu}\varphi(\eta)\tag{103}$$

the resultant equation along with the boundary conditions can be solved with different numerical approaches such as Runge-Kutta, Runge-Kutta-Fehlberg, compact finite difference, high-order finite difference; however, in this present paper, the Thomas algorithm with finite difference discretization will be used to solve the Hiemenz profile in the Julia framework.

#### 2.2.2. Numerical Solution of Hiemenz Flow Problem

The computational approach for the Hiemenz flow is similar to the Blasius solution since the final ordinary equation of the Hiemenz flow (Equation (98)) is similar to the Blasius similarity solution (Equation (20)). If it is assumed that *ϕ*- = *h* and *ϕ* = *p* in Equation (98), the system of equations will be:

$$h'' + ph' - h^2 + 1 = 0\tag{104}$$

$$p' - h = 0,\tag{105}$$

where the boundary conditions of the system are:

$$p(0) = 0\tag{106}$$

$$h(0) = 0\tag{107}$$

$$h(\infty) = 1.\tag{108}$$

If finite difference scheme is applied to the *h* variable as in Blasius flow. The final system of equations will be:

$$A\_n h\_{n+1} + B\_n h\_n + C\_n h\_{n-1} + D\_n = 0,\tag{109}$$

where the *An*, *Bn*, *Cn*, and *Dn* are:

$$A\_{\rm nl} = \frac{1}{\Delta \eta^2} + \frac{p\_{\rm n}}{2 \Delta \eta} \tag{110}$$

$$B\_{\rm ll} = \frac{-2}{\Delta \eta^2} - h\_{\rm ll} \tag{111}$$

$$C\_n = \frac{1}{\Delta \eta^2} - \frac{p\_n}{2\Delta \eta} \tag{112}$$

$$D\_n = 1.\tag{113}$$

The relation between *hn* and *hn*−<sup>1</sup> is taken as it is done in Blasius flow. When Equation (37) is substituted into the Equation (109), the final *Gn* and *Hn* coefficients can be found as:

$$G\_{\rm II} = \frac{-\mathbb{C}\_{\rm n} G\_{n-1} + D\_{\rm n}}{B\_{\rm n} + \mathbb{C}\_{\rm n} H\_{\rm n-1}} \tag{114}$$

$$H\_{\rm II} = \frac{-A\_{\rm II}}{B\_{\rm II} + \mathbb{C}\_{\rm n} H\_{\rm n-1}}.\tag{115}$$

In order to start the calculation, *H*<sup>1</sup> can be assumed as 0 and the *G*<sup>1</sup> can be calculated as 0 from the boundary conditions. Once *G*<sup>1</sup> and *H*<sup>1</sup> are calculated, *Gn* and *Hn* can be calculated from Equations (114) and (115). *p* can be calculated as it is done for Blasius flow:

$$p\_n = p\_{n-1} + \int\_{\eta\_{n-1}}^{\eta\_n} h d\eta. \tag{116}$$

The system of equations can be implemented in the Julia environment as it is shown in Listing 4. The linear profile assumption is taken as the initial condition of *h* (see Figure 2) as it is done for the Blasius solution. The reason *hp* and *errorPro file* variables are used in the Listing 4 is that the linear profile converges to the Hiemenz profile in each iteration so in order to check the difference between previous and present profiles, the solution vector from the previous iteration is copied and it is compared with the new solution vector. If the difference between these two solution vectors is less than , which is an arbitrary limit and can be taken as 1 <sup>×</sup> <sup>10</sup>−8, then it can be said that the solution has converged. It has to be noted that can be taken as any number; however, if it is small, the solution will be more accurate. The final result of the Hiemenz flow can be seen in Figure 7. The results are also validated with White [44].

#### 2.2.3. Homann Flow Problem

Homann flow is similar to Hiemenz flow. The only difference is that Homann flow is an axisymmetric version of the Hiemenz flow. The same schematic (Figure 6) can represent this flow as well; however, in this flow, the circle is the projection of a sphere. On the other hand, it was the projection of a cylinder in Hiemenz flow. Derivation of the Homann similarity solution has the almost same procedure as well but it uses cylindrical coordinates instead of Cartesian coordinates. The velocity components of the flow can be shown as:

$$
\upsilon\_r = \upsilon\_r(r, z) \tag{117}
$$

$$w\_{\theta} = 0\tag{118}$$

$$
\upsilon\_z = \upsilon\_z(r, z) \tag{119}
$$

$$p = p(r, z). \tag{120}$$

**Figure 7.** The velocity distribution of Hiemenz similarity solution obtained by Julia code and data digitized from White [44].

As axisymmetric assumptions, derivative with respect to *θ* ( *<sup>∂</sup> ∂θ* ) and *v<sup>θ</sup>* can be assumed as zero. The Navier–Stokes equations in cylindrical coordinates can be written as:

$$\frac{\partial v\_{\ell}}{\partial r} + \frac{v\_{r}}{r} + \frac{\partial v\_{z}}{\partial z} = 0 \tag{121}$$

$$v\_{r}\frac{\partial v\_{r}}{\partial r} + v\_{z}\frac{\partial v\_{z}}{\partial z} = -\frac{1}{\rho}\frac{\partial p}{\partial r} + \nu \left(\frac{\partial^{2}v\_{r}}{\partial r^{2}} + \frac{1}{r}\frac{\partial v\_{r}}{\partial r} - \frac{v\_{r}}{r^{2}} + \frac{\partial^{2}v\_{r}}{\partial z^{2}}\right) \tag{122}$$

$$
v\_r \frac{\partial v\_z}{\partial r} + v\_z \frac{\partial v\_z}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial z} + \nu \left( \frac{\partial^2 v\_z}{\partial r^2} + \frac{1}{r} \frac{\partial v\_z}{\partial r} + \frac{\partial^2 v\_z}{\partial z^2} \right). \tag{123}
$$

If the same procedures applied to Navier–Stokes in Cartesian coordinates for Hiemenz flow, are applied in cylindrical coordinates, the final potential function and stream function will be:

$$
\phi(r, z) = k(\frac{r^2}{2} - z^2) \tag{124}
$$

$$
\psi(r, z) = kzr^2,\tag{125}
$$

where *k* is a constant which depends on the freestream flow and the body shape. The corresponding velocities are:

$$
\omega\_r = \frac{\partial \phi}{\partial r} = \frac{1}{r} \frac{\partial \psi}{\partial z} = kr \tag{126}
$$

$$
v\_z = \frac{\partial \phi}{\partial z} = \frac{-1}{r} \frac{\partial \psi}{\partial r} = -2kz. \tag{127}$$

Once velocities are calculated, the pressure from the Bernoulli can be calculated as:

$$p\_0 - p = \frac{1}{2}\rho k^2 (r^2 + 4z^2). \tag{128}$$

The viscous velocity components and pressure can be calculated as:

$$w\_{\mathbb{T}} = r \mathbb{g}(z) \tag{129}$$

$$v\_z = -2f(z)\tag{130}$$

$$p\_0 - p = \frac{1}{2}\rho k^2 (r^2 + F(z)).\tag{131}$$

**Listing 4.** Implementation of Thomas algorithm for Hiemenz profile. ✞ ☎

```
1 h = [(i−1)/(N−1) for i=1:N]
2
3 for i = 2:N
4 p[i] = p[i−1] + (h[i] + h[i−1])*Δη/2
5 end
6
7 while 1e−8<=errorProfile
8
9 A = [ 1/Δη2 + p[i]/(2*Δη) for i=1:N]
10 B=[−2/Δη2 − h[i] for i=1:N]
11 C = [ 1/Δη2 − p[i]/(2*Δη) for i=1:N]
12 D=[1 for i=1:N]
13
14 for i=2:N−1
15 G[i] = − ( C[i]*G[i−1] + D[i] )/(B[i] + C[i] * H[i−1])
16 H[i] = − A[i] /(B[i] + C[i] * H[i−1])
17 end
18
19 hp = copy(h)
20
21 for i=N−1:−1:2
22 h[i] = G[i] + H[i] * h[i+1]
23 end
24
25 errorProfile = maximum( abs .(hp−h))
26
27 for i = 2:N
28 p[i] = p[i−1] + (h[i] + h[i−1])*Δη/2
29 end
30 end
```
✝ ✆ After the affine transformation, the final ordinary differential equation will be:

$$
\varphi^{\prime\prime\prime} + 2\varphi\varphi^{\prime\prime} - \varphi^{\prime 2} + 1 = 0,\tag{132}
$$

where the boundary conditions of the equation is:

$$
\varphi'(0) = 0\tag{133}
$$

$$
\varphi(0) = 0\tag{134}
$$

$$
\varphi'(\infty) = 1\tag{135}
$$

If one compares Hiemenz and Homann similarity solutions (see Equations (98) and (132)) the only difference is the 2 in the second term. In the same manner, the computational process will be the same except for two lines of code.

#### 2.2.4. Numerical Solution of Homann Flow Problem

Since equations of Hiemenz and Homann similarity solutions are the same except one coefficient, procedures for the solution are also the same, except for two lines of code. If it is assumed that *ϕ*-= *h* and *ϕ* = *p* in Equation (132), the system of equations will be:

$$h'' + 2ph' - h^2 + 1 = 0\tag{136}$$

$$p'-h=0,\tag{137}$$

where the boundary conditions of the system are:

$$p(0) = 0\tag{138}$$

$$h(0) = 0\tag{139}$$

$$h(\infty) = 1\tag{140}$$

If the finite difference scheme is applied to the *h* variable, the final system of equations will be:

$$A\_n h\_{n+1} + B\_n h\_n + C\_n h\_{n-1} + D\_n = 0,\tag{141}$$

where the *An*, *Bn*, *Cn*, and *Dn* are:

$$A\_{n} = \frac{1}{\Delta \eta^{2}} + \frac{p\_{n}}{\Delta \eta} \tag{142}$$

$$B\_n = \frac{-2}{\Delta \eta^2} - h\_n \tag{143}$$

$$C\_n = \frac{1}{\Delta \eta^2} - \frac{p\_n}{\Delta \eta} \tag{144}$$

$$D\_{\text{ll}} = 1.\tag{145}$$

The final *Gn* and *Hn* coefficients are the same as Hiemenz flow as well and they are:

$$G\_{\rm n} = \frac{-\mathbb{C}\_{\rm n} G\_{\rm n-1} + D\_{\rm n}}{B\_{\rm n} + \mathbb{C}\_{\rm n} H\_{\rm n-1}} \tag{146}$$

$$H\_{\rm II} = \frac{-A\_{\rm II}}{B\_{\rm II} + \mathbb{C}\_{\rm n} H\_{\rm n-1}}.\tag{147}$$

The final code can be seen in Listing 5. If the Hiemenz code (4) and Homann code (5) are compared, the only difference is in the *A* and *C* and it is because of the finite difference approach for the *h*- .

**Listing 5.** Implementation of Thomas algorithm for Homann profile. ✞ ☎

```
1 h = [(i−1)/(N−1) for i=1:N]
2
3 for i = 2:N
4 p[i] = p[i−1] + (h[i] + h[i−1])*Δη/2
5 end
6
7 while 1e−8<=errorProfile
8
9 A = [ 1/Δη2 + p[i]/Δη for i=1:N]
10 B=[−2/Δη2 − h[i] for i=1:N]
11 C = [ 1/Δη2 − p[i]/Δη for i=1:N]
12 D=[1 for i=1:N]
13
14 for i=2:N-1
15 G[i] = − ( C[i]*G[i−1] + D[i] )/(B[i] + C[i] * H[i−1])
16 H[i] = − A[i] /(B[i] + C[i] * H[i−1])
17 end
18
19 hp = copy(h)
20
21 for i=N−1:−1:2
22 h[i] = G[i] + H[i] * h[i+1]
23 end
24
25 errorProfile = maximum( abs .(hp−h))
26
27 for i = 2:N
28 p[i] = p[i−1] + (h[i] + h[i−1])*Δη/2
29 end
30 end
```
The final solution profile can be seen in Figure 8. The results are also validated with White [44]. One can use these results to validate their own codes. Reference data will be shared on GitHub as well.

**Figure 8.** The velocity distribution of Homann similarity solution obtained by Julia code and digitized from White [44].

#### **3. Discussion**

Fast computational fluid dynamics solvers are crucial for engineers because the design process requires a lot of simulations to reach the final and optimized design. One of the most important factors that make a solver fast is the language itself. Writing almost the same script in different languages may give different solution central processing unit (CPU) times. For instance, a script that uses long and complex for loops in both Julia and Python environments will result in different solution times because Python is slower with for loops in general. The reasons that cause this slow behavior will not be covered in this paper since it is out of the scope of this paper. However, it is important to state these differences in order to decide which coding language is proper for the simulation that is planned.

Sometimes, some other criteria, such as a user-friendly environment, might be the critical condition. The easy matrix–matrix multiplication and backslash linear system solver can be some of the user-friendly examples. In CFD, the usage of vectors and matrices is so common. Most of the time, it is required to multiply or add vectors or matrices with one another. Fortran, which is one of the fastest languages and also one of the fundamental languages in the CFD industry, does not have a built-in element-wise vector or matrix multiplication feature. This requires the use of for loops for each element-wise matrix and vector operation. As a result of this, every time, one needs to use a for loop to do these operations. This ends up with hard-to-read codes and also excessive usage of indices is another source of possible mistakes that will cause trouble during the debugging. However, in the Julia environment, it can be done in one line without any trouble. Listing 6 is showing the two for loop usage. Both of them use Julia syntax in order to prevent confusion; however, it is required to state that Fortran has a different syntax than this but, logically, the for loops are the same with Fortran logic.

Another user-friendly feature of Julia over Fortran is to plot a vector in the code with built-in functions. However, in Fortran, it is not possible to do it with built-in functions so one needs to extract the solution vector or matrix to an external file to visualize it. As it is stated before, coding language selection can be an important topic. In this section, the strong sides of the Julia environment will be stated. One of the user-friendly features of Julia is compact for loop syntax. Listing 7 shows the traditional way and compact way to use a for loop in the Julia environment.

**Listing 6.** Comparison of the element-wise and compact matrix-matrix multiplication. ✞ ☎

```
1 for l=1:m
2 for k=1:nz
3 for j=1:ny
4 for i=1:nx
5 A[i,j,k,l] = B[i,j,k,l]*C[i,j,k,l]
6 end
7 end
8 end
9 end
10
11 A .= B.*C
```
**Listing 7.** Comparison of the traditional and compact way of for loop usage. ✞ ☎

```
1 for i=1:N
2 h[i] = (i−1)/(N−1)
3 end
4
5 h = [(i−1)/(N−1) for i=1:N] ✝ ✆
```
The readability of the code is also extremely important in order to explain the code to others or provide it as open-source code. Julia is so strong in this topic because it allows the usage of LATEX language within the code. It is one of the unique abilities of this language. If you are writing a formula that is full of Greek letters and regular letters, it is inevitable to use complex variable names. Comparison of the usage of Greek letters and regular formula writing can be seen in Listing 8.

✝ ✆

**Listing 8.** Comparison of the Greek letter usage and conventional formula writing in Julia. ✞ ☎

```
1 alpha = a*beta−c/(gamma−1)*(d_i−1/(eta^2))
2
3 α = a*β−c/(γ−1)*(di−1/(η2)) ✝ ✆
```
The last advantage of Julia stated in this present paper is the ability to do both dynamic and static RAM allocations. It may be important to manually allocate the variable sizes and types for optimized code. However, sometimes it is easier to use just variables without any initializing. Initializing differences in the Julia environment can be seen in Listing 9.

**Listing 9.** Different initializing methods in Julia environment. ✞ ☎

```
1 nx = 10
2 ny = 15.2
3
4 q = Array{Float64}(undef , nx)
5 for i=1:nx
6 q[i] = i
7 end
8
9 a = Array{Float64}(undef , nx)
10 c = zeros (nx)
11 a=q
12 b=q
13 c=q
✝ ✆
```
As it is seen from the example scripts, Julia is a user-friendly, fast, open-source, and free language which can increase productivity drastically [36]. It is a great choice both for those who are new to coding and coding experts. Julia also can call the C, Fortran, and Python libraries so it is great for experienced engineers who think their previous code in other coding languages will be useless.

#### **4. Conclusions**

In the computational fluid dynamics industry, it is crucial to have some predictions about the flow that will be simulated. It helps to spot the location at which finer mesh is required. Fundamental knowledge about canonical flows is crucial in this point because most

of the complex flow consists of a combination of a couple of canonical flows. For example, an airfoil CFD simulation can consist of boundary-layer flow over the smooth part of it, mixing layer flow where the tail ends, and blunt body flow in the wake region. In other words, one airfoil simulation consists of three different canonical flows. A full understanding on canonical flows is extremely important in order to simulate similar flows accurately and cheaply. In this paper, boundary-layer theory is introduced and boundary-layer flows are derived from scratch. The Blasius flow, Hiemenz flow, Homann flow, Falkner–Skan flow are the focus of this paper. Once the derivations of them are completed, derived forms are implemented in the Julia environment. In order to model the equations, a finite difference scheme for space discretization is used and the Thomas algorithm is used for the linear system solution. It has to be noted that other methods such as Runge-Kutta, Runge-Kutta-Fehlberg, compact finite difference, high-order finite difference can also be used to solve the given ODEs; however, in this present paper, authors preferred the finite difference method and Thomas algorithm. The authors further discussed the advantages of the Julia language over some other coding languages available in the literature. It is shown that Julia syntax is so straightforward, easy, and user-friendly. Strong sides of the Julia environment are stated with the given comparisons as well. The popularity of Julia may drastically increase in the near future because of its potential.

**Author Contributions:** code generation, F.O. and K.K.; validation, F.O. and K.K.; writing—original draft preparation, F.O.; writing—review and editing, F.O. and K.K.; visualization, F.O. All authors have read and agreed to the published version of the manuscript.

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

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All the data used and generated in this study is available in the GitHub link provided in Appendix A.

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

#### **Appendix A**

Julia setup files can be downloaded from their website (https://julialang.org/downloads/ (accessed on 1 June 2021)). The website also includes instructions on how to install Julia on Windows, Linux, and mac operating systems. Some of the useful resources for learning Julia are listed below:


It is common to use external packages for Julia. In order to do that, Pkg, which is Julia's built-in package manager, can be used. Once Julia is opened, Pkg can be activated with the "]" button in Windows. In Linux, calling "julia" in the terminal will open it. After that "Pkg.add("Pluto")" will trigger the setup process for that package. In here, we used Pluto as an example because, in GitHub, our codes are developed in the Pluto environment. After Pluto is installed. Pluto can be run with "Pluto.run()". This command will open a new tab in the browser which you can run your Julia codes. After that, the "using Pluto" line must be placed to the top of the file. For "Plots" package, the commands will be "Pkg.add("Plots")" and "using Plots". Since the Plots package does not have a GUI, there is not a command called "Plots.run()".

Other than Pluto, JuliaPro which includes Julia and the Juno IDE (https://juliacomputing. com/products/juliapro/ (accessed on 1 June 2021)) can be used as an editor and compiler. This software contains a set of packages for plotting, optimization, machine learning, database, and much more. Pluto is appropriate for small scripts while JuliaPro is better for more complex codes. The GitHub link of the codes used in this paper is:

• https://github.com/frkanz/A-CFD-Tutorial-in-Julia (accessed on 1 June 2021)

#### **References**


## *Article* **PyDA: A Hands-On Introduction to Dynamical Data Assimilation with Python**

#### **Shady E. Ahmed, Suraj Pawar and Omer San \***

School of Mechanical and Aerospace Engineering, Oklahoma State University, Stillwater, OK 74078, USA; shady.ahmed@okstate.edu (S.E.A.); supawar@okstate.edu (S.P.)

**\*** Correspondence: osan@okstate.edu; Tel.: +1-405-744-2457; Fax: +1-405-744-7873

Received: 1 November 2020; Accepted: 26 November 2020; Published: 29 November 2020

**Abstract:** Dynamic data assimilation offers a suite of algorithms that merge measurement data with numerical simulations to predict accurate state trajectories. Meteorological centers rely heavily on data assimilation to achieve trustworthy weather forecast. With the advance in measurement systems, as well as the reduction in sensor prices, data assimilation (DA) techniques are applicable to various fields, other than meteorology. However, beginners usually face hardships digesting the core ideas from the available sophisticated resources requiring a steep learning curve. In this tutorial, we lay out the mathematical principles behind DA with easy-to-follow Python module implementations so that this group of newcomers can quickly feel the essence of DA algorithms. We explore a series of common variational, and sequential techniques, and highlight major differences and potential extensions. We demonstrate the presented approaches using an array of fluid flow applications with varying levels of complexity.

**Keywords:** data assimilation; variational and sequential methods; Kalman filtering; forward sensitivity; measurements fusion

#### **1. Introduction**

Data assimilation (DA) refers to a class of techniques that lie at the interface between computational sciences and real measurements, and aim at fusing information from both sides to provide better estimates of the system's state. One of the very mature applications that significantly utilize DA is weather forecast, that we rely on in our daily life. In order to predict the weather (or the state of any system) in the future, a model has to be solved, most often by numerical simulations. However, a few problems rise at this point and we refer to only few of them here. First, for accurate predictions, these simulations need to be initiated from the true initial condition, which is never known exactly. For large scale systems, it is almost impossible to experimentally measure the full state of the system at a given time. For example, imagine simulating the atmospheric or oceanic flow, then you need to measure the velocity, temperature, density, etc. at every location corresponding to your numerical grid! Even in the hypothetical case when this is possible, measurements are always contaminated by noise, reducing the fidelity of your estimation. Second, the mathematical model that completely describes all the underlying processes and dynamics of the system is either unknown or hard to deal with. Then, approximate and simplified models are adopted instead. Third, the computational resources always constrain the level of accuracy in the employed schemes and enforce numerical approximations. Luckily, DA appears at the intersection of all these efforts and introduces a variety of approaches to mitigate these problem, or at least reduce their effects. In particular, DA techniques combines possibly incomplete dynamical models, prior information about initial system's state and parameterization, and sparse and corrupted measurement data to yield optimized trajectory in order to describe the system's dynamics and evolution.

As indicated above, dynamical data assimilation techniques have a long history in computational meteorology and geophysical fluid dynamics sciences [1–3]. Then comes the question: why do we write such an introductory tutorial about a historical topic? Before answering this question, we highlight a few points. DA borrows ideas from numerical modeling and analysis, linear algebra, optimization, and control. Although these topics are taught separately in almost every engineering discipline, their combination is rarely presented. We believe that incorporating DA course in engineering curricula is important nowadays as it provides a variety of global tools and ideas that can potentially be applied in many areas, not just meteorology. This was proven while administering a graduate class on "Data Assimilation in Science and Engineering" at Oklahoma State University, as students from different disciplines and backgrounds were astonished by the feasibility and utility of DA techniques to solve numerous inverse problems they are working on, not related to weather forecast. Nonetheless, the availability of beginner-friendly resources has been the major shortage that students suffered from. The majority of textbooks either derives DA algorithms from their very deep roots or surveys their historical developments, without focus on actual implementations. On the other hand, available packages are presented in a sophisticated way that optimizes data storage and handling, computational cost, and convergence. However, this level of sophistication takes a steep learning curve to understand the computational pipeline as well as the algorithmic steps, and a lot of learners fall hopeless during this journey.

Therefore, the main objective of this tutorial paper is to familiarize beginning researchers and practitioners with basic DA ideas along with easy-to-follow pieces of codes to feel the essence of DA and trigger the priceless "aha" moments. With this in mind, we choose Python as the coding language, being a popular, interpreted language, and easy to understand even with minimum programming background, although not the most computationally favored language in high performance computing (HPC) environments. Moreover, whenever possible, we utilize the built-in functions and libraries to minimize the user coding efforts. In other words, the provided codes are presented for demonstrative purposes only using an array of academic test problems, and significant modifications should be incorporated before dealing with complex applications. Meanwhile, this tutorial will give the reader a jump-start that hopefully shortens the learning curve of more advanced packages. A Python-based DA testing suite has been also designed to compare different methodologies [4].

Dynamical data assimilation techniques can be generally classified into variational DA and sequential DA. Variational data assimilation which works by setting an optimization problem defined by a cost functional along with constraints that collectively incorporate our knowledge about the system. The minimizer of the cost functional represents the DA estimate of the unknown system's variables and/or parameters. On the other hand, in sequential methods (also known as statistical methods), the system state is evolved in time using background information until observations become available. At this instant, an update (correction) to the system's variables and/or parameters is estimated and the solver is re-initialized with this new updated information until new measurements are collected, and so on. We give an overview of both approaches as well as basic implementation. In particular, we briefly discuss the three dimensional variational data assimilation (3DVAR) [5,6], the four dimensional variational data assimilation (4DVAR) [7–12], and forward sensitivity method (FSM) [13,14] as examples of variational approaches. Kalman filtering and its variants [15–22] are the most popular applications of sequential methods. We introduce the main ideas behind standard Kalman filter and its extensions for nonlinear and high-dimensional problems. The famous Lorenz 63 is utilized to illustrate the merit of all presented algorithms, being a simple low-order dynamical system that exhibit interesting dynamics. This is to help readers to digest the different pieces of codes and follow the computational pipeline. Then, the paper is concluded with a section that provides the deployment of selected DA approaches for dynamical systems with increasing levels of dimensionality and complexity. We highlight here that the primary purpose of this paper is to provide an introductory tutorial on the data assimilation for educational purposes. All Python

implementations of the presented algorithms as well as the test cases are made publicly accessible at our GitHub repository https://github.com/Shady-Ahmed/PyDA.

#### **2. Preliminaries**

#### *2.1. Notation*

Before we dive into the technical details of dynamical data assimilation approaches, we briefly present and describe our notations and assumptions. In general, we assume that all vector-valued functions or variables are written as a column vector. Unless stated otherwise, boldfaced lowercase letters are used to denote vectors and boldface uppercase letters are reserved to matrices. We suppose that the system state at any time *<sup>t</sup>* is denoted as **<sup>u</sup>**(*t*)=[*u*1(*t*), *<sup>u</sup>*2(*t*), ... , *un*(*t*)]*<sup>T</sup>* <sup>∈</sup> <sup>R</sup>*n*, where *<sup>n</sup>* is the state-space dimension. The dynamics of the system are governed by the following differential equation

$$\frac{d\mathbf{u}}{dt} = \mathbf{f}(\mathbf{u}; \boldsymbol{\theta}),\tag{1}$$

where **<sup>f</sup>** : <sup>R</sup>*<sup>n</sup>* <sup>×</sup> <sup>R</sup>*<sup>p</sup>* <sup>→</sup> <sup>R</sup>*<sup>n</sup>* encapsulates the model's dynamics, with *<sup>θ</sup>* <sup>∈</sup> <sup>R</sup>*<sup>p</sup>* being the vector of model's parameters and *p* being the number of these parameters. With a time-integration scheme applied, the discrete-time model can be written as follows,

$$\mathbf{u}(t\_{k+1}) = M(\mathbf{u}(t\_k); \boldsymbol{\theta}), \tag{2}$$

where *M* is the one-time step transition map that evolves the state at time *tk* to time *tk*<sup>+</sup><sup>1</sup> = *tk* + Δ*t*, with Δ*t* being the time step length.

We denote the true value of the state variable as **u***t*, which is assumed to be unknown and a good approximation of it is sought. Our prior information about the state **u** is called the background, with a subscript of *b* as **u***b*. This represents our beginning knowledge, which might come from historical data, numerical simulations, or just an intelligent guess. The discrepancy between this background information and true state is denoted as *ξ<sup>b</sup>* = **u***<sup>t</sup>* − **u***b*, resulting from imperfect model, inaccurate model's initialization, incorrect parameterization, numerical approximations, etc. From probabilistic point of view, we suppose that the background error has a zero mean and a covariance matrix of **<sup>B</sup>**. This can be represented as *<sup>E</sup>*[*ξb*] = 0 and *<sup>E</sup>*[*ξbξ<sup>T</sup> <sup>b</sup>* ] = **<sup>B</sup>**, where **<sup>B</sup>** <sup>∈</sup> <sup>R</sup>*n*×*<sup>n</sup>* is a symmetric and positive-definite matrix and the superscript *T* refers to the transpose operation. Moreover, we assume that unknown true state has a multivariate Gaussian distribution with a mean **u***<sup>b</sup>* and a covariance matrix **B** (i.e., **u***<sup>t</sup>* = N (**u***b*, **B**)).

We define the set of the collected measurements at a specific time *tk* as **<sup>w</sup>**(*tk*) <sup>∈</sup> <sup>R</sup>*m*, where *<sup>m</sup>* is the dimension of observation-space. We highlight that the observed quantity need not be the same as the state variable. For instance, if the state variable that we are trying to resolve is the temperature of sea surface, we may have access only to radiance measurements by satellites. However, those case be related to each other through Planch–Stefan's law, for instance. Formally, we can relate the observables and the state variables as

$$\mathbf{w}(t\_k) = h(\mathbf{u}\_t(t\_k)) + \mathbf{\mathcal{f}}\_{m'} \tag{3}$$

where *<sup>h</sup>* : <sup>R</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup>*<sup>m</sup>* defines the mapping from state-space to measurement-space and *<sup>ξ</sup><sup>m</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* denotes the measurement noise. The mapping *h* can refer to the sampling (and probably interpolation) of state variables at the measurements locations, relating different quantities of interest (e.g., relating see surface temperature to emitted radiance), or both! The model's map *M* and observation operator *h* can be linear, nonlinear, or a combination of them. Similar to the background error, the observation noise *ξ<sup>m</sup>* is assumed to possess a multivariate normal distribution, with a zero mean and a covariance matrix **<sup>R</sup>** <sup>∈</sup> <sup>R</sup>*m*×*m*, i.e., *<sup>ξ</sup><sup>m</sup>* <sup>=</sup> <sup>N</sup> (0, **<sup>R</sup>**). An extra grounding assumption is that the measurement noise and the state variables (either true or background) are uncorrelated. Furthermore, all noises are assumed to be temporally uncorrelated (i.e., white noise). Even though we consider only Gaussian distribution for the background error, and observation noise, we emphasize that there has been a lot of studies dealing with non-Gaussian data assimilation [23–26].

The objective of data assimilation is to provide an algorithm that fuses our prior information **u***<sup>b</sup>* and measurement data **w** to yield a better approximation of the unknown true state. This better approximation is called the analysis, and denoted as **u***a*. The difference between this better approximation and the true state is denoted as *ξ<sup>a</sup>* = **u***<sup>t</sup>* − **u***a*.

#### *2.2. Twin Experiment Framework*

In a realistic situation, the true state values are unknown and noisy measurements are collected by sensing devices. However, for testing ideas, the ground truth need to be known beforehand such that the convergence and accuracy of the developed algorithm can be evaluated. In this sense, the concept of twin experiment has been popular in data assimilation (and inverse problems, in general) studies. First, a prototypical test case (all called toy problems!) is selected based on the similarities between its dynamics and real situations. Similar to your first "Hello World!" program, the Lorenz 63 and Lorenz 96 are often used in numerical weather forecast investigations, the one-dimensional Burgers equation is explored in computational fluid dynamics developments, the two-dimensional Kraichnan turbulence and three-dimensional Taylor–Green vortex are analyzed in turbulence studies, and so on. A reference true trajectory is computed by fixing all parameters and running the forward solver until some final time is reached. Synthetic measurements are then collected by sampling the true trajectory at some points in space and time. A mapping can be applied on the true state variables and arbitrary random noise is artificially added (e.g., a white Gaussian noise). Finally, the data assimilation technique of interest is implemented starting from false values of the state variables or the model's parameters along with the synthetic measurement data. The output trajectory of the algorithm is thus compared against the reference solution, and the performance can be evaluated. It is always recommended that researchers get familiar with twin experiment frameworks as they provide well-structured and controlled environments for testing ideas. For instance, the influence of different measurement sparsity and/or level of noise can be cheaply assessed, without the need to locate or modify sensors.

#### **3. Three Dimensional Variational Data Assimilation**

The three dimensional variational data assimilation (3DVAR) framework can be derived from either an optimal control or Bayesian analysis points of view. The interested readers can be referred to other resources for mathematical foundations (e.g., [27]). In order to compute a good approximation of the system state, the following cost functional can be defined,

$$f(\mathbf{u}) = \frac{1}{2}(\mathbf{w} - h(\mathbf{u}))^T \mathbf{R}^{-1} (\mathbf{w} - h(\mathbf{u})) + \frac{1}{2}(\mathbf{u} - \mathbf{u}\_b)^T \mathbf{B}^{-1} (\mathbf{u} - \mathbf{u}\_b),\tag{4}$$

where the first term penalizes the discrepancy between the actual measurement **w** and the state variable mapped into the observation space *h*(**u**) (also called the model predicted measurement). The second term aims at incorporating the prior information, weighted by the inverse of the covariance matrix to reflect our confidence in this background. We highlight that all terms in Equation (4) are evaluated at the same time, and thus the 3DVAR can be referred to as a stationary case.

The minimizer of *J*(**u**) (i.e., the analysis) can be obtained by setting the gradient of the cost functional to zero as follows,

$$\nabla f(\mathbf{u}) = -\mathbf{D}\_h^T(\mathbf{u}\_d)\mathbf{R}^{-1}(\mathbf{w} - h(\mathbf{u}\_d)) + \mathbf{B}^{-1}(\mathbf{u}\_d - \mathbf{u}\_b) = 0,\tag{5}$$

where **<sup>D</sup>***h*(**u**) <sup>∈</sup> <sup>R</sup>*m*×*<sup>n</sup>* is the Jacobian matrix of the operator *<sup>h</sup>*(**u**). The difficulty of solving Equation (5) depends on the form of *h*(**ua**) as it can either be linear, or highly nonlinear.

#### *3.1. Linear Case*

For linear observation operator (i.e., *h*(**u**) = **Hu**, and **D***h*(**u**) = **H**, where **H** is an *m* × *n* matrix), the evaluation of the analysis **u***<sup>a</sup>* in Equation (5) reduces to solving the following linear system of equations

$$(\mathbf{B}^{-1} + \mathbf{H}^T \mathbf{R}^{-1} \mathbf{H})\mathbf{u}\_a = (\mathbf{B}^{-1} \mathbf{u}\_b + \mathbf{H}^T \mathbf{R}^{-1} \mathbf{w}).\tag{6}$$

We note that (**B**−<sup>1</sup> <sup>+</sup> **<sup>H</sup>***T***R**−1**H**) on the left-hand side is an *<sup>n</sup>* <sup>×</sup> *<sup>n</sup>* matrix, and hence this is called the model-space approach to 3DVAR. Furthermore, a popular incremental form can be derived from Equation (6) by adding and subtracting **H***T***R**−1**Hu***<sup>a</sup>* to/from the right-hand side and rearranging to get the following form,

$$\mathbf{u}\_a = \mathbf{u}\_b + (\mathbf{B}^{-1} + \mathbf{H}^T \mathbf{R}^{-1} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{R}^{-1} (\mathbf{w} - \mathbf{H} \mathbf{u}\_b). \tag{7}$$

Moreover, the Sherman–Morrison–Woodbury inversion formula can be used to derive an observation-space solution to the 3DVAR problem (for details, see [27], page 327) as follows,

$$\mathbf{u}\_{a} = \mathbf{u}\_{b} + \mathbf{B} \mathbf{H}^{T} (\mathbf{R} + \mathbf{H} \mathbf{B} \mathbf{H}^{T})^{-1} (\mathbf{w} - \mathbf{H} \mathbf{u}\_{b}). \tag{8}$$

Note that (**<sup>R</sup>** <sup>+</sup> **HBH***T*) is an *<sup>m</sup>* <sup>×</sup> *<sup>m</sup>* matrix, compared to (**B**−<sup>1</sup> <sup>+</sup> **<sup>H</sup>***T***R**−1**H**) being an *<sup>n</sup>* <sup>×</sup> *<sup>n</sup>* matrix. Thus, Equation (7) or Equation (8) might be computationally favored based on the values of *n* and *m*. We also highlight that in either cases, matrix inversion is rarely (almost never) computed directly, and efficient linear system solvers should be utilized, instead. An example of a Python function for the implementation of the 3DVAR algorithm with a linear operator is shown in Listing 1.

**Listing 1.** Implementation of 3DVAR for with a linear observation operator. ✞ ☎

```
import numpy as np
def Lin3dvar(ub,w,H,R,B,opt):