Suite-CFD: An Array of Fluid Solvers Written in MATLAB and Python

Abstract

Computational Fluid Dynamics (CFD) models are being rapidly integrated into applications across all sciences and engineering. CFD harnesses the power of computers to solve the equations of fluid dynamics, which otherwise cannot be solved analytically except for very particular cases. Numerical solutions can be interpreted through traditional quantitative techniques as well as visually through qualitative snapshots of the flow data. As pictures are worth a thousand words, in many cases such visualizations are invaluable for understanding the fluid system. Unfortunately, vast mathematical knowledge is required to develop one’s own CFD software and commercial software options are expensive and thereby may be inaccessible to many potential practitioners. To that extent, CFD materials specifically designed for undergraduate education are limited. Here we provide an open-source repository, which contains numerous popular fluid solvers in $2D$ (projection, spectral, and Lattice Boltzmann), with full implementations in both MATLAB and Python3. All output data is saved in the $.vtk$ format, which can be visualized (and analyzed) with open-source visualization tools, such as VisIt or ParaView. Beyond the code, we also provide teaching resources, such as tutorials, flow snapshots, measurements, videos, and slides to streamline use of the software.

1. Introduction

Computational Fluid Dynamics (CFD) models are being applied to problems across all sciences and engineering. From designing more aerodynamic sportswear and vehicles, to understanding animal locomotion, to personalized medicine, to disease transmission, and to predicting hurricanes and atmospheric phenomena, it is difficult to find situations where greater knowledge of the underlying fluid dynamics isn’t desired or valuable. Due to its immense importance, many scientists have dedicated their careers to the field, whether they study particular fluid phenomena or develop tools, either experimental or numerical, for other scientists and engineers to use. Possibly because of how vital understanding fluid dynamics is to human flourishing, proving existence, uniqueness, and the possible breakdown of solutions to the system of non-linear partial differential equations that govern fluid dynamics is one of the Millennium Problems, to which the Clay Institute offers a 1 million dollar prize for solving [1]. On the other hand, CFD allows us to solve these equations, the Navier-Stokes equations, which are the equations that detail the conversation of momentum and mass for a fluid. For an incompressible, viscous fluid, they can be written as follows,

$$\rho \left[\frac{\partial \mathbf{u}}{\partial t}(\mathbf{x},t)+\left(\mathbf{u}(\mathbf{x},t)·\nabla \right)\mathbf{u}(\mathbf{x},t)\right]=-\nabla p(\mathbf{x},t)+\mu \Delta \mathbf{u}(\mathbf{x},t)$$

(1)

$$\nabla ·\mathbf{u}(\mathbf{x},t)=0$$

(2)

where $\mathbf{u}(\mathbf{x},t)$ and $p(\mathbf{x},t)$ are the fluid’s velocity and pressure, respectively. The physical properties of the fluid are given by $\rho$ and $\mu$, its density and dynamic viscosity, respectively. The independent variables are the time, t, and spatial position, $\mathbf{x}$. Here all variables that pertain to the fluid, e.g., $\mathbf{u}$ and p, are written in an Eulerian framework on a fixed Cartesian mesh, $\mathbf{x}$. Think of a fixed Cartesian mesh as a rectangular grid to which over the course of the simulation physical quantities, like velocity, pressure, etc., are measured at specific lattice points. In an Eulerian framework, specific blobs of fluid are not tracked over time, as though they are billiard balls, but rather, we focus on measuring quantities of interest at specific locations in space, located on a fixed grid. Note that Equations (1) and (2) provide the conservation of momentum and mass, respectively.

Many techniques have been developed to solve the above system of equations, such as projection methods [2,3,4,5,6], spectral methods [7,8,9], and Lattice Boltzmann methods [10,11,12]; however, they are typically left out of undergraduate curricula due to the mathematical background required for implementation, while they are commonly found in graduate curricula [13]. Moreover, it is not only mathematical knowledge that is a barrier; CFD requires expertise in at least four knowledge domains [14]: flow physics, numerical methods, computer programming, and validation, either experimental or theoretical. Thankfully, this has not stopped academics from developing and researching methods to integrate CFD into undergraduate courses [13,14,15,16,17,18,19,20,21]. Some have offered best practices guides and advice for future CFD undergraduate courses [13,21]. Stern et al. [13] has suggested for beginners to focus on the overall CFD process and flow visualizations, in order to solidify their fundamental understanding of fluid physics. Appropriate flow visualizations also help students validate and verify their simulations through careful inspection [21].

As computer literacy is rapidly becoming more of a central focus in undergraduate curricula, CFD could be a natural outlet where students may not only implement and test algorithms, but also practice proper data analysis and data visualization, while building computational experience. Much of the fluid dynamics software that has been used in previous courses is either foreign to students or only commercially available [14,17,19,20]. Commercial software has been recommended by the American Physical Society and American Association of Physics Teachers as a tool to help students prepare for 21st century careers by exposing them to a greater range of experiences and opportunities [22,23]. Unfortunately due to expensive software licenses, some schools, particularly smaller institutions, may not be able to afford such software, especially if they would only be used intermittently to complement the existing course material.

Recently, popular programming languages selected for many undergraduate curricula in science, mathematics, and engineering have been either MATLAB or Python [24,25,26,27,28]. MATLAB [29] and Python [30] are both interpreted languages, making it easier (and more attractive) for students to begin programming [31]. Moreover if students have already been exposed to these programming languages, giving them a platform to refine their skills in these languages while at the same strengthening their intuition in fluid dynamics may prove beneficial.

For this reason, many scientists, engineers, and mathematicians who are passionate about education have begun developing CFD software intended for educational purposes written in these languages [32,33,34,35,36]. Such open source software packages can be integrated into courses in a well-streamlined fashion, with tutorials, activities, exercises, and notes provided, and may even blur the line between classroom activities and contemporary research [32,35,36,37]. In particular L.A. Barba and G. Forsyth [33] have developed a guide for students to numerically solve the Navier-Stokes equations in 12 scaffolding steps and L.A. Barba and O. Mesnard [34] have developed a similar scaffolding lesson structure for students to study classical aerodynamics using potential flow models. Both are written in Python using Jupyter Notebooks [38]. Pawar and San [39] developed modules for teaching advanced undergraduate and graduate students how to develop their own fluid solvers in the Julia programming language [40].

Rather than develop a collection of modules that teach how to implement particular CFD numerical schemes, we offer the scientific community a variety of popular fluid solvers that solve traditional problems in fluid dynamics, with two independent but equal implementations written in MATLAB and Python. The emphasis is not placed on students having to develop or implement numerical schemes, but to allow them the opportunity to get an accelerated start in CFD, and run, tweak, and analyze simulations of traditionally studied problems in fluids courses, as well as explore data visualization and how to present data that is appropriate for fluid physics interpretation. All of this can be done while testing their own hypotheses in familiar programming environments. In the remainder of this manuscript, we will give a high-level overview of the CFD schemes currently implemented in the software (Section 2), guided tutorials on how to run, visualize, and analyze simulations (Section 3), and provide numerous examples to which students may elect to study or modify as well (Section 4).

Furthermore, we provide an in-depth mathematical description of each numerical scheme in Appendix B. The open-source software repository can be accessed at https://github.com/nickabattista/Holy_Grail. Any simulation data produced from the software can be visualized and analyzed in open-source software, such as VisIt [41], maintained by Lawrence Livermore National Laboratory, or ParaView [42], developed by Kitware, Inc., Los Alamos National Laboratory, and Sandia National Labs. We also provide teaching resources, such as slides, figures, and movies in the Supplementary Materials.

2. Brief Overview of the Three Fluid Solvers

In this section we will give a brief overview of the three fluid solvers currently implemented in the software—a projection method, a spectral (FFT) method, and the Lattice-Boltzmann method. Our goal is to provide students a high-level overview of the methods and how they compare to one another. Further details regarding their mathematical foundations and implementations are given in Appendix B.

2.1. Projection Method

The projection method was first introduced by Chorin in 1967 [2] and independently by Temam in 1968 [4], to solve the viscous, incompressible Navier-Stokes equations. Projection methods are finite difference based numerical schemes [43], in which the fluid equations are solved in an Eulerian framework, i.e., the computational grid is static and fluid dynamics variables, like velocity, pressure, etc., are measured at particular locations on the grid over time during a simulation. Projection methods have been extensively used throughout the fluid dynamics community for decades, in numerous applications across many fields, while being continually improved for efficiency and accuracy [5,44,45,46].

In a nutshell, this method decouples the velocity and pressure fields, using operator splitting and a Helmholtz-Hodge decomposition. This makes it possible to explicitly solve Equations (1) and (2) in a few steps [2] while also increasing computational efficiency. Please see Appendix B.1 for further mathematical details.

Furthermore, projection methods allow one to define explicit boundary conditions (BCs) for velocity on edges of the computational domain, whether they are Dirichlet, Neumann, or Robin boundary conditions [5,47]. In our software, we allow users to impose tangential velocity boundary conditions along the edges of the rectangular domain, while the normal direction boundary conditions for velocity are explicitly set to zero. The choice to set normal BCs to zero was made in order to simplify the code by ensuring that volume conservation would be automatically satisfied, i.e., situations will not arise in which a specified amount of fluid is being pumped into the domain while a different amount is leaving the domain per unit time.

Having defined the fluid velocity as $\mathbf{u}=(u,v)$, one can set the boundary conditions for u on the top and bottom of the domain ($v=0$ along these edges) or v on the left and right sides ($u=0$ along these edges). See Figure 1 for an illustration. Later in Section 4, we will showcase two different traditional problems in fluid dynamics with a projection method:

• Cavity Flow
• Circular Flow in a Square Domain

The examples are different due to the changes in boundary conditions that the user can impose, see Figure 1. In Figure 1 we illustrate the boundary conditions for each example: cavity flow (a) and circular flow (b). This figure also illustrates that the domain could be constructed to be either rectangular or square and that the boundary conditions do not have to be uniform among all sides, nor across the domain.

We will also compare different fluid scales in both examples, using the Reynolds Number, $Re$. The Reynolds Number is given by

$$Re=\frac{\rho LV}{\mu },$$

(3)

where $\rho$ and $\mu$ are the fluid’s density and dynamic viscosity, respectively, while L and V are a characteristic length- and velocity-scale of the system. In essence, $Re$ is a ratio of inertial forces to viscous forces in a system. If $Re>>1$, inertia dominates, if $Re<<1$, viscous forces dominate, and if $Re=1$, the inertia and viscous forces are balanced. An example of high $Re$ flow would be the flow around a marble quickly falling through air. An example of low $Re$ flow would be the fluid flow around a marble slowly falling through honey.

2.2. Spectral Method (FFT)

Spectral methods are a different class of differential equation solvers. They are well-known for achieving highly accurate solutions, i.e., spectral accuracy [48,49,50]. Spectral accuracy occurs when error decreases exponentially with a small change in grid resolution, e.g., there is a linear relationship between the logarithm of error and grid resolution, i.e.,

$$log\left(error\right)\sim N,$$

where N is the number of grid points. These methods approximate solutions using orthogonal expansions, such as Fourier Series or other orthogonal basis of functions, like Chebyshev Polynomials or Legendre Polynomials. In contrast to finite difference based schemes, like the projection method, where solutions are approximated at specific locations in space, in spectral methods a particular orthogonal series expansion is chosen and the earnest is on finding the coefficients of the expansion.

Here we chose to use the basis of Discrete Fourier Transform (DFT). The DFT can be used to express functions that are not periodic, unlike Fourier Series expansions, which are used to represent periodic functions. Moreover, we used the Fast Fourier Transform (FFT), which yields the same numerical values as a DFT, but is considerable more computationally efficient, i.e., fast.

Thus, we implemented a FFT based spectral scheme to solve the viscous, incompressible Navier-Stokes equations in $2D$. We also chose to solve these equations in their vorticity formulation (see Appendix B.2 for mathematical details). The vorticity formulation will lead to a few numerical benefits.

In a nutshell, this will recast our problem into solving for the vorticity, $\mathit{\omega }=\omega \widehat{k}$, and the streamfunction, $\psi$. Velocity can be defined in terms of the curl of the streamfunction, i.e., a vector potential, see below

$$\mathbf{u}=\nabla \times \psi \widehat{k}.$$

(4)

Hence this allows us to find the velocity field, $\mathbf{u}=(u,v)$, once we find $\psi$, e.g.,

$$u=\frac{\partial \psi }{\partial y}\mathrm{and}v=-\frac{\partial \psi }{\partial x}.$$

Moreover, the incompressibility condition (Equation (2)) will be automatically satisfied since $\mathbf{u}$ is defined to be the curl of a scalar, and the divergence of the curl of a vector is always identically zero.

Furthermore, working in the vorticity formulation requires us to only solve a Poisson equation for the streamfunction in terms of the vorticity, $\mathit{\omega }$, i.e.,

$$\Delta \psi =-\omega .$$

(5)

Notice that Equation (5) is a linear equation. Moreover, upon taking the FFT of the Equation (5), the computation is transformed into frequency (Fourier) space, which offers two advantages—increased speed and accuracy.

In practice, we will solve for the streamfunction at the next time-step based on the previously computed vorticity and then update the vorticity using information from the newly updated streamfunction. We chose to use a Crank-Nicholson time-stepping scheme [51] to update the vorticity to the next time-step. The Crank-Nicholson scheme is well-known for being unconditionally-stable for diffusion problems and as well as for being 2nd order accurate in time and space, as an implicit method [52]. These accuracy and stability properties make the Crank-Nicholson an ideal candidate for numerically solving such equations.

Although using a FFT (spectral method) preserves high accuracy of the spatial information at a particular time-step, numerical error still unfortunately creeps in due to the time-stepping nature of solving an evolution equation for the fluid’s momentum (see Equation (A15) in Appendix B.2). That is, evolving the system forward in time is where the brunt of the error takes place, rather than the spatial discretization.

In contrary to the projection method, this FFT-based approach allows us to explore periodic boundary conditions. For example, if fluid flow is moving vertically-upwards through the top of the computational domain, it will re-enter moving vertically-upwards through the bottom of the domain. This has advantages and disadvantages. Some advantages are that explicit boundary conditions do not need to be satisfied nor enforced, which aids in computational speed-up. However, to that extent, one disadvantage is that this makes it more difficult to enforce particular boundary conditions when desired. FFT-based fluid solvers have been used in the fluid-structure interaction community, in particular within immersed boundary methods [32,37,53], when one wants to study the interactions of an object and the fluid to which it’s immersed. As the focus is on the object and fluid interactions, it is typically modeled away from the domain boundaries in the middle of the computational domain to avoid boundary artifacts. For a mathematically detailed description of this FFT-based spectral method see Appendix B.2.

For this particular FFT-based fluid solver for the vorticity-formulation of the viscous, incompressible Navier-Stokes equations, we will highlight the following examples:

• Side-by-Side Vorticity-Region Interactions
• Evolution of Vorticity from an Initial Velocity Field

Contrary to the examples shown for the projection method in Section 2.1, these examples are not differentiated by different boundary conditions, but rather different initial vorticity configurations in the computational domain. Figure 2 gives the initial vorticity configurations for the two cases considered: side-by-side vorticity region interactions (left) and an initial vorticity field that defined by a velocity vector field (right). In Figure 2a, counterclockwise (CCW) and clockwise (CW) correspond to regions of uniform vorticity, where vorticity initialized as a positive or negative constant for $CCW$ and $CW$, respectively. In all other regions of the computational domain, the vorticity is either initialized to zero (Figure 2a). Although not shown, one could also initialize simulations on a rectangular grid. In Figure 2b, the initial vorticity is defined by computing the vorticity from a simulation’s velocity vector field. The velocity vector field was taken from a time-point in a Pulsing_Heart simulation [36] in the open-source IB2d software [32,37]. In this example, we illustrate that if a snapshot of the velocity field is known, it is possible to evolve the fluid dynamics forward using this spectral (FFT) method, even though it is based on the vorticity formulation of the Navier-Stokes equations.

2.3. Lattice Boltzmann Method

The Lattice Boltzmann method (LBM) does not explicitly (or implicitly) solve the incompressible, Navier-Stokes equations, rather it uses discrete Boltzmann equations to model the flow of a fluid [11]. In a nutshell, the LBM tracks fictitious particles of fluid flow, thinking of the problem more as a transport equation, e.g.,

$$\frac{\partial f}{\partial t}+\mathbf{u}·\nabla f=\Omega ,$$

(6)

where $f(\mathbf{x},t)$ is the particle distribution function, i.e., a probability density, u is the fluid particle’s velocity, and $\Omega$ is what is called the collision operator. However, rather than these particles moving in a Lagrangian framework, the Lattice Boltzmann method simplifies this assumption and restricts the particles to the nodes on a lattice.

In traditional CFD methods, one typically discretizes the fluid equations in an Eulerian framework, as in a projection or spectral method. When you are free of such restriction, as in the LBM, it allows one to more readily deal with complex boundaries, model porous structures and multi-phase flow, as well as incorporate microscopic interactions, while also allowing for massive parallelization of the algorithm, leading to increased computational efficiency [12,54,55,56].

The LBM involves collision and streaming steps, where fictional fluid particles consecutively propagate (and collide) over a discretized lattice, and the fluid density is evolved forward. There are multiple ways to handle boundary conditions in the LBM. One such way is to use bounce-back boundary conditions, which are executed by masking points in the domain via boolean expressions, thus defining regions where the fluid can and cannot flow [55]. In short, the propagating (streaming) directions are simply reversed when they hit a boundary node, see Figure A3 in Appendix B.3. The bounce-back conditions can be used to enforce the necessary and physical no-slip boundary conditions for fluid flow as well as defining solid objects in the interior of the grid (for an example see Figure 33). For a detailed mathematical description of the LBM see Appendix B.3.

To highlight the LBM fluid solver, we showcase the following examples:

• Flow past one or more cylinders
• Flow past porous cylinder

3. How to Run the Simulations, Visualize, and Analyze

In this section we will give instructions on how to run, visualize, and analyze the simulations produced in this software. We will explain how to change parameters below as well. In particular, we will provide both overviews of what is possible to visualize and analyze as well as a guided walk-through of the steps required, designated by Guide in their title. We broke this section into the following subsections:

• Running a Simulation
• Visualizing the Data in VisIt
• Guide: Running the Spectral (FFT) Method’s ‘bubble3’ Example
• Guide: Visualizing the Spectral (FFT) Method’s ‘bubble3’ Data
• Guide: Analyzing the Spectral (FFT) Method’s ‘bubble3’ Data

3.1. Running a Simulation

Here we will briefly describe how to run each fluid solver’s corresponding simulations. To illustrate the software structure and how to run the simulations, we will use the MATLAB (https://www.mathworks.com/products/matlab.html) [29] version. Note that the Python implementation is consistent in its structure as well as naming conventions.

   To run the simulation, one would need to do the following:

• Either clone the Holy_Grail repository or download the Holy_Grail zip file at https://github.com/nickabattista/Holy_Grail/ to your local machine. Note you can download or clone this repository to any directory on your local machine.
• Open MATLAB and go to the appropriate sub-directory for the particular fluid solver example you wish to run, e.g.,
  • Projection_Method: for the projection method solver and its examples
  • FFT_NS_Solver: for the spectral (FFT) method solver and its examples
  • Lets_Do_LBM: for the Lattice Boltzmann method solver and its examples
  and go into the MATLAB subfolder. (If you wanted to run an example in Python, you would change directories until you are in the complementary Python subfolder.)
• The main script name for each fluid solver is named accordingly, e.g.,
  • Projection_Method.m - main script for the projection method
  • FFT_NS_Solver.m - main script for the spectral (FFT) solver
  • lets_do_LBM.m - main script for the Lattice Boltzmann solver
  Each of these scripts contains all the functions and modules necessary to run each example. The file print_vtk_files.m produces the $.vtk$ data files during each simulation. The user should not have to change this file, unless they do not wish to save some of the data for storage reasons. If this is the case, this can be accomplished by commenting out the lines corresponding to whichever data is not to be saved.
• Depending on the subfolder that you are in, to run an example, you would type its main script name into the MATLAB command window and click enter.
• The simulations are not instantaneous; they may take a few minutes. Upon starting the simulation, a folder named vtk_data is produced, which will be used for storing the simulation data files in $.vtk$ format. As the simulation runs, it will dump more data into this folder. Note that depending on the simulation you choose and the grid resolution, the simulations may take from a few minutes to ∼30 min. Moreover, if you wish to run multiple simulations, note that the default folder name is set to vtk_data, and so running additional simulations will overwrite any previously saved data in the previous vtk_data folder. Please manually change the name of the folder after each simulation to avoid this.

Digging deeper into each fluid solver script:

• If one desires to change parameters for a particular simulation, e.g., grid resolution, viscosity (either $\mu$ (Projection) or $\nu$ (FFT) or $\tau$ (LBM), etc.), options are available to the user without digging deep into any of the sub-functions or individual modules beyond the main function itself, see Figure 3.
• Figure 3 highlights the main choices the user can make for each of the three fluid solvers (a) Projection, (b) Spectral (FFT), and (c) Lattice Boltzmann. The user can identify which of the built-in examples they wish to run by changing the appropriate string flag labeled choice.
• Note that each example will run as described in Section 4 if all the parameters are selected to be those in the parameter Table that provided within each example’s section. Caution: numerical stability conditions may be violated by choice of other parameters. To test other parameter values, we suggest making incremental variations, e.g., try doubling or halving viscosity, rather than testing values orders of magnitude apart. Numerical instability will appear if solvers cannot converge or begin dumping NaNs (not-a-number), which represent unrepresentable numbers to computers, due to their inherent floating-point arithmetic.
• As the entire fluid-solver is contained within each script, the user can go deeper into the code to modify examples or create their own. This can be done in the following sub-functions:
  • Projection Method: please_Give_Me_BCs (choice) 
  • Spectral (FFT) Solver: please_Give_Initial_Vorticity_State (choice, NX, NY) 
  • Lattice Boltzmann Solver: give_Me_Problem_Geometry (choice, Nx, Ny, percentPorosity)

3.2. Visualizing the Data in VisIt

In this section we will briefly describe how produce visualizations, such as those Section 4. We will briefly detail the steps to visualize the data using the open-source software VisIt (https://visit.llnl.gov) [41]. Note that each simulation produces data in the $.vtk$ format and so one could alternatively use the open-source software ParaView [42] for visualization purposes as well. Later, in Section 3.3 we will guide a user through running an example (the spectral (FFT) method’s bubble3 example) followed by a step-by-step guide to visualizing its corresponding data in Section 3.4. We used VisIt version 2.13.3.

To visualize the data, one would need to do the following:

• Open VisIt [41]
• Open the desired Eulerian data (magnitude of velocity, vorticity, velocity vector field, etc.) from the vtk_data folder.
• A table of the possible data stored and their corresponding name when saved is found below in 

  Table 1. Data stored to $.vtk$ format and its corresponding storage name.

  Parameter	Type	‘.vtk’ Name
  Magnitude of Velocity	Scalar	$uMag$
  Velocity Vectors	Vector	u
  Horizontal Velocity	Scalar	$uX$
  Vertical Velocity	Scalar	$uY$
  Pressure (projection only)	Scalar	P
  Vorticity	Scalar	$Omega$
  Geometry (LBM only)	Mesh	$Bounds$

  .
• To Visualize Geometry (LBM only):

  (a)

  Click Open

  (b)

  Go to the vtk_data data folder that the simulation produced

  (c)

  Click on the grouping of Bounds, click OK

  (d)

  In VisIt, click on Add then Mesh→mesh.

  (e)

  Then click Draw

  (f)

  You can elect to change the color of boundary or size by double clicking on the Mesh in the VisIt database listing window.

• To Visualize Velocity Vectors:

  (a)

  Click Open

  (b)

  Go to the vtk_data data folder that the simulation produced

  (c)

  Click on the grouping of u, click OK

  (d)

  In VisIt, click on Add then Vector→u

  (e)

  Then click Draw. An error message might pop up during at time-step zero (first data point) if velocity is identically zero

  (f)

  You can elect to change the color, size, scaling, and number of velocity vectors by double-clicking on u in the VisIt database listing window

  (g)

  To add streamlines of the velocity field, first make sure that your Active Source is set to u (see Figure 4a), then click Add→Pseudocolor→operators→IntegralCurve→u (see Figure 4b). Note you have the option to put numerous seed points for streamlines to stem in the domain. Here we have chosen to use a line of points, sampled it at 4 evenly spaced locations within that line, and specified particular integration tolerances, and a maximum number of steps (see Figure 4c). Click Draw. You can also adjust the streamline aesthetics to how you desire, e.g., color, thickness, etc.

• To Visualize the Eulerian scalar data (e.g., Vorticity, Magnitude of Velocity, etc.):

  (a)

  Click Open

  (b)

  Go to the vtk_data data folder that the simulation produced

  (c)

  Click on the grouping of the desired Eulerian scalar data, for example, Omega (for Vorticity), click OK

  (d)

  In VisIt, click on Add then Pseudocolor→Omega.

  (e)

  Then click Draw

  (f)

  You can elect to change the colormap and/or colormap scaling by double clicking on Omega in the VisIt database listing window. The following table (

  Table 2. Colormaps used for each scalar variable for visualization.

  Parameter	Colormap
  Magnitude of Velocity	RdYlBu or BrBG (Section 4.1 only)
  Horizontal Velocity	RdYlBu
  Vertical Velocity	RdYlBu
  Pressure	orangehot (Section 4.1 only)
  Vorticity	hot_desaturated
  FTLE	PuRd

  ) details the specific colormap used for each scalar variable in the manuscript,

  (g)

  To visualize the finite-time Lyapunov exponent (FTLE), your Active Source in Visit (see Figure 5a) must be on the velocity vectors, u. Then click Add→Pseudocolor→operators→LCS→u (see Figure 5b). We used the attributes in Figure 5c for computing the FTLE. In particular, we limited the maximum advection time to 0.05 and maximum number of steps to 10. Once those are changed, click Draw.

  (h)

  To add contours of the desired scalar variable, first make sure that your Active Source (see Figure 5a) is set to the correct variable, then click Add→Contour→<desired variable name>. Note that you have the option to scale the contours separate from the colormap of the variable itself. Moreover, for FTLE your active source must be u and you can modify how FTLE are computed, e.g., Figure 5c; however, we used consistent values for both FTLE computations.

3.3. Guide: Running the Spectral (FFT) Method’s ‘bubble3’ Example

In this guided walk-through, we wish to illustrate how to run a simulation and change its parameters while also observing what information is being printed to the screen. We will use the Spectral (FFT) Method’s bubble3 example and run two different examples where each uses a different fluid kinematic viscosity, $\nu$.

To run the simulation(s), do the following:

• First we must make sure that we are inside of MATLAB the corresponding directory for the MATLAB version of the spectral-based solver. To do this click through: Holy_Grail→FFT_NS_Solver→MATLAB, see Figure 6. Note that I have downloaded (or cloned) the software into my Documents folder.
• To run this simulation, type FFT_NS_Solver into MATLAB’s command window and click enter. See Figure 7 for what should print to the command window shortly thereafter. Note that the bubble3 example is the default built-in example upon downloading the software (see Figure 3b).
  Upon starting the simulation, information regarding the simulation solver and simulation is printed to the screen. Moreover, the current_time within the simulation gets printed to the screen at the time-points in which the data is stored. The data is stored in the vtk_data folder, which is made upon starting the simulation.
• Once the simulation has finished running, we will change the vtk_folder’s name to Simulation_1_Data. This folder contains all of the data that was stored during the simulation, i.e., the vorticity, fluid velocity field, magnitude of velocity, etc. It is imperative to change this folder name after each simulation if more simulations are to be run, or else the new data will printed over the previous simulation’s data.
• Once the folder’s name has been changed, we will open the FFT_NS_Solver.m script to change one (or more) of the parameters. Here we will only change the fluid’s kinematic viscosity, $\nu$; however, there are other options in which you could change, such as the computational domain size and/or grid resolution, see Figure 3b. As an example change $\mathtt{nu}$ from $1.0\times {10}^{-3}\to 1.0\times {10}^{-2}$, see Figure 8. Once you have changed $\nu$, you can close the script.
  Note that if you chose to change the grid resolution and domain size, make sure that the resolution in x and y are equal, i.e., $dx=L_x/N_x=L_y/N_y=dy$. Moreover if you increase the resolution (increase $N_x,N_y$) the simulations will take longer to run and the resulting data will require more storage.
• Finally, run this new simulation case by typing FFT_NS_Solver into the command window and pressing enter. Once that simulation has finished running, change its folder name to Simulation_2_Data. Now that we have a couple of simulations run, and have their corresponding data saved, we will next focus on visualizing each simulation’s dynamics using the open-source software VisIt [41].

3.4. Guide: Visualizing the Spectral (FFT) Method’s ‘bubble3’ Data

Here we will guide the user through visualizing the bubble3 simulation data from Section 3.3. In particular we will visualize the magnitude of velocity data. Note that other Eulerian scalar data, like vorticity, the x- or y-component of velocity, etc., may be visualized analogously, as we are following the protocol from the previous section on data visualization, Section 3.2’s Visualizing the Eulerian scalar data. To visualize the fluid’s velocity vector field, please follow Section 3.2’s steps for To Visualize Velocity Vectors.

To visualize the uMag data, do the following:

• First we must open the desired magnitude of velocity data. Recall that the fluid’s magnitude of velocity data was saved under the name uMag. Click open in the VisIt toolbar and in the dialog box that pops up, locate the desired data directory, here it is Simulation_Data_1 from earlier, and select the uMag group. Figure 9 is provided as a visual aid for these steps.
• Although we have opened the data, nothing will appear visualized, yet. However, now uMag is the Active Source in the window, see Figure 10. We must now tell VisIt how to visualize the data. We will visualize magnitude of velocity as a background colormap, as it is a scalar quantity. To do this click ‘Add’ and then select Pseudocolor from the menu that appears, followed by uMag, i.e., Add→Pseudocolor→uMag.
• To finally visualize the data, click Draw. Window 1 will then show the fluid’s magnitude of velocity data at time-point 0, i.e., the initial magnitude of velocity, see Figure 11.
• Furthermore, we note that certain choices of colormaps are better for both conveying the data, but also for people who are colorblind. Generally, the default colormap, which was used in Figure 11, is a bad choice, as it is impossible for people who are red-green colorblind to interpret the field [57]. Thus, we will change the colormap to another, e.g., the RdYlBu (red-yellow-blue) colormap, by double-clicking on the uMag.*.vtk database:Pseudocolor uMag in the database list and appropriately selecting the desired colormap. Other properties may also be changed in this box, such as the colormap’s scaling.
  Also, you can visualize the magnitude of velocity at the different time-points that were saved during the simulation by moving the time-slider. Figure 12 shows the data at time-point 30. Moreover, the colormap and its properties may be modified by double-clicking on uMag.*.vtk database:Pseudocolor uMag in the database list.
• Note that the time in the simulation does not generally correspond to the time-point stored. For example, in this case the simulation modeled 30 s of time and 60 total data time-points were stored. Thus, data was saved every 0.5 s of simulation time. Furthermore, you could repeat Steps 1–4 for data in Simulation_2_Data from the second simulation you ran. Instead of comparing the 30th time-point’s data, we will compare the 60th, i.e., the last time-point. Repeating this entire process for the other simulation’s data yields the comparison as shown in Figure 13, below.
  However, we cannot qualitatively compare the magnitude of velocity between each of these simulations as they are in Figure 13. Each simulation’s colormap has a different scale. In order to compare, we must make both scales equivalent.
• To change the scale, double-click on each of the databases. This will bring up a dialog window (as described earlier) where you can change the scaling, e.g., minimum and maximum, as well as the colormap itself. Change the minimum and maximum values to 0.0 and 0.2, respectively. Figure 14 illustrates where to change the minimum and maximum for the colormap’s scale. Once those values are changed, click Apply. Note that you must do this for both of the databases individually.
• After changing the scale, the comparison is significantly different than when first visualized, see Figure 15.
  Furthermore, we also provide visualizations of other time-points and/or scalar data in Appendix C.

3.5. Guide: Analyzing the Spectral (FFT) Method’s ‘bubble3’ Data

To analyze the data in VisIt, we will provide a walk-through of the steps that necessary to measure a fluid quantity (particular variable) across a line of interest within the computational domain using the Spectral (FFT) Method’s bubble3 example from Section 3.3. We will assume that you have visualized the data for the magnitude of velocity that corresponds to both of the simulations. In this tutorial we will compare the magnitude of velocity for both simulations across a vertical line through the center of the computational domain at the last time-point stored, i.e., the 60th time-point. See Figure 16 for an idea of where the measurement will take place.

To analyze the uMag data, do the following:

• Make sure that Active source is the uMag data corresponding to the data from Simulation_1_Data in VisIt, see Figure 16.
  Once the Active Source has been appropriately selected, in the command bar, in the VisIt menu go to Controls→Query (see Figure 16).
• Next select Lineout under Queries and then under Variables select Scalars and then uMag, i.e., Variables→Scalars→uMag, see Figure 17.
• Next we will define the line segment in which we intend to measure the data across. To do this the user gets to choose the starting point and ending point of the line segment, see Figure 17b. We will measure across a vertical line through the center of the domain, hence we choose $(x,y,z)=(0.5,0.0,0.0)$ and $(x,y,z)=(0.5,1.0,0.0)$ as the starting and ending point, respectively. We also will measure the value of uMag at 250 sampled points. Thus, change Sample Points to 250, see Figure 17b and select Use sampling. VisIt will automatically interpolate the data so that it is measured at number of equally-spaced sample points that is designated by the user.
• Once you click Query a plot should pop up in a new Active Window in Visit, i.e., Active window 2, see Figure 18. If you double-click on the Lineout(uMag) bar in the VisIt database listing window, you can change various attributes about the line, such as its thickness or color among others. Moreover, you can use the time-slider to observe how the magnitude of velocity across the chosen line varies from time-point to time-points.
• To get back to the colormap window, select Active Window 1, rather than 2.
• You can repeat this procedure and overlay multiple uMag measured data on the same plot, each corresponding to a different simulation, e.g., repeat this process for Simulation_2_Data. Be sure to choose the correct data for the Active source.
• If you repeat this entire process for the Simulation_2_Data and move the time-slider such that the data being shown in Active window 2 corresponds to the last time-point stored (the 60th), you should see what is give in Figure 19. Also, note that we have changed each line’s thickness individually.

4. Built-in Examples (for Each Fluid Solver)

In this section we will highlight a subset of the built-in examples that can be run immediately upon download (or git clone) of the software. These examples were selected as they are either traditional problems in fluid dynamics or problems that could naturally lead to fruitful discussion and analysis of the underlying dynamics. The examples we will discuss are:

• Lid-Driven Cavity Flow via the Projection Method
• Circular Flow in a Square Domain via the Projection Method
• Side-by-side Vorticity Region Interactions via the Spectral (FFT) Method
• Evolution of Vorticity from an Initial Velocity Field via the Spectral (FFT) Method
• Flow Past One or More Cylinders via the Lattice Boltzmann Method
• Flow Past a Porous Cylinders via the Lattice Boltzmann Method

For each of these examples, we will provide the necessary simulation details, including a brief background about problem as well as the numerical parameters that were used in each simulation. Furthermore, we will also make observations regarding each simulation’s dynamics, investigate some of the data, and suggest future ideas for how students could either modify or explore each example further.

These examples are available within both the MATLAB and Python implementations in the software. Note that Section 3 described how to run, visualize, and analyze the simulations, and thus the visualizations and analysis contained here can be recreated by students and other practitioners.

4.1. Cavity Flow (via Projection Method)

In this example we will use the Projection Method in the software to run a cavity flow problem. Note that this example is contained in the Projection Method folder and may be run by selecting the ‘cavity_top’ option (see line 58 in Figure 3a).

For this lid-driven cavity flow problem, the non-zero horizontal velocity on the top wall, $u=u_T$, is set to $u_T=4.0$ m/s, while all other tangential (and normal) velocities along the boundaries are set to zero (see Figure 1a). Note that the top wall velocity is not immediately set at $u_T$, but rather the flow ramps up from 0 to the preset $u_T$ along this wall, to avoid instantaneous acceleration artifacts. For the simulations shown below with $Re=4000$, use the computational parameters listed in

Table 3. Numerical parameters for the Cavity simulation for $Re=4000$.

Parameter	Variable	Units	Value
Domain Size	$[L_x,L_y]$	m	$[1,2]$
Spatial Grid Resolution	$[N_x,N_y]$		$[128,256]$
Spatial Grid Size	$dx=dy$	m	$L_x/N_x=L_y/N_y$
Time Step Size	$dt$	s	${10}^{-3}$
Total Simulation Time	T	s	6
Fluid Density	$\rho$	$\mathrm{kg}/{\mathrm{m}}^3$	1000
Fluid Dynamic Viscosity	$\mu$	$\mathrm{kg}/(\mathrm{m}·\mathrm{s})$	1

. The Reynolds Number was set by defining the characteristic length scale to be $L=L_x$, the horizontal length of the domain, and the characteristic velocity to be $V=u_T$, e.g., $Re=\frac{\rho L_x{u}_T}{\mu }$. To change the Reynolds Number, the dynamic viscosity was varied. For $Re=4000,400,40$ and 4, we used $\mu =1,10,100$, and $1000\mathrm{kg}/(\mathrm{m}·\mathrm{s})$ , respectively. Note that for the cases with $Re=4$ and 40, we decreased the time-step to $dt=5\times {10}^{-5}\mathrm{s}$ to ensure numerical stability of solutions. If we had not, errors would have been magnified every time-step of the simulation and/or the numerical solutions may have blown up. In CFD it is common that one may need to vary the time-step depending on the $Re$ or other parameters of the system. Situations in which very small time-steps need to be taken to ensure numerical stability of solutions are known as stiff equations, for more information please see [43].

Upon having run the cavity model for $Re=4000$, we visualized its corresponding simulation data using VisIt [41], as illustrated by Figure 20, which gives colormaps of vorticity, magnitude of velocity, velocity in the horizontal direction throughout the domain, and pressure, as well as depictions of the velocity field, both as a scaled and non-scaled quantity. For the scaled velocity field snapshot, the scale factor was equal to the largest magnitude of velocity across the entire computational domain. In the non-scaled plot, streamlines are also given. Streamlines show the path that a passive particle would take in the flow at a particular moment in time. The data shown here was taken at $t=6.0\mathrm{s}$, when the simulation ended.

A fully-formed vortex developed by that time near the top of the cavity, while a smaller vortex appears to be forming below (see Vorticity pane in Figure 20). Snapshots illustrating the formation of such vortices from this case of $Re=4000$ are presented in Figure 21, which gives a colormap of vorticity and the velocity vector field. As fluid is moving along the top edge of the domain from left-to-right, the fluid nearest to the top begins moving in the same direction. Eventually fluid down in the cavity begins moving in the same direction as well, until it reaches the right boundary, where it must move downwards. It is easier for the fluid to move downward because of the faster flows moving towards the right above it, due to the boundary condition. As the fluid moves downward along the wall, fluid towards the middle of the cavity begins moving downwards as well. In tandem, with the fluid moving left-to-right along the top and top-to-bottom along the right, these fluid patterns initiate the formation of a clockwise-spinning vortex.

Now that there is a vortex spinning clockwise, it causes the fluid below the vortex (moving right-to-left) to begin moving right-to-left until it reaches the wall on the left and an analogous situation occurs to the above, except an oppositely-spinning (counter-clockwise) vortex begins to form. Recall that these vortices can both be traced back to the fluid moving along the top of the domain left-to-right. Within the region that the second vortex forms, much of the energy (velocity) has dissipated away, resulting in a smaller, less strong vortex.

We can explicitly measure the horizontal velocity descending down the cavity. This data is presented in Figure 22, which gives the horizontal velocity vs. depth in the cavity during four different snapshots in the simulation. We also measure the horizontal velocity across three different lines descending into the cavity. Near the surface of the cavity (at depth = $0\mathrm{m}$), the velocity is equal to the boundary condition. As one descends into the cavity, there are spikes in both the positive and negative horizontal velocities. These spikes correspond to locations near the edges of the vortices, where there is faster moving fluid.

Finally, one can change the $Re$ by varying the dynamic viscosity, $\mu$ to observe how vortex formation changes within the cavity, as well as, how quickly energy gets dissipated in higher viscosity settings. Recall that if $\mu$ is changed to 40 or 4 the time-step, $dt$, was adjusted for as described above. Figure 23 illustrates how vortex formation changes within the cavity for different $Re$ at three different snapshots during the simulations. Moreover, it gives the horizontal velocity vs. cavity depth, as measured across the middle of the cavity. Note that the $Re=4$ and $Re=40$ data virtually overlaps in the figure, and velocity (and system’s energy) quickly dissipates to zero going further into the cavity.

Students may elect to try the following:

• Recreate the above results for cavity flow with the geometric and simulation parameters given in Table 3 for one or more $Re$.
• Change the domain width (cavity width) to observe differences as the cavity gets wider for a given $Re$.
• Vary the lid-velocity to observe differences (e.g., varying $Re$ for different velocities as opposed to differing viscosities, $\mu$).

4.2. Circular Flow in a Square Domain (via Projection Method)

In this example we will use the Projection Method in the software to run an example involving circular flow in a square domain. Note that this example is in the Projection Method folder and may be run by selecting the ‘whirlwind’ option (see line 58 in Figure 3a).

In this case, all the tangential boundary conditions along the computational domain were set to be $U=u_{top}=-u_{bot}=v_{left}=-v_{right}=1.0$ m/s, see Figure 1b. As mentioned in Section 4.1 above, the tangential wall velocities are not immediately set to U, but rather the flow ramps up along each edge from 0 to the preset $u_{top},u_{bot},v_{left}$, or $v_{right}$ along this wall, to avoid instantaneous acceleration artifacts. For the simulation shown below with $Re=4000$, use the computational parameters listed in

Table 4. Numerical parameters for the whirlwind simulation for $Re=4000$.

Parameter	Variable	Units	Value
Domain Size	$[L_x,L_y]$	m	$[1,1]$
Spatial Grid Resolution	$[N_x,N_y]$		$[256,256]$
Spatial Grid Size	$dx=dy$	m	$L_x/N_x=L_y/N_y$
Time Step Size	$dt$	s	${10}^{-3}$
Total Simulation Time	T	s	24
Fluid Density	$\rho$	$\mathrm{kg}/{\mathrm{m}}^3$	1000
Fluid Dynamic Viscosity	$\mu$	$\mathrm{kg}/(\mathrm{m}·\mathrm{s})$	$0.25$

. Note that as $Re=\frac{\rho LU}{\mu }$, with $L=L_x=L_y$ and $U=1.0$ m/s, one can vary the dynamic viscosity, $\mu$, to change the $Re$. If you decide to make $Re$ smaller (e.g., increasing $\mu$), you may also need to decrease the time-step, $dt$, significantly to ensure numerical stability, as mentioned in Section 4.1. The requirement that very small time-steps are necessary to ensure numerical stability denotes what are called stiff equations. For more information on stiff equations, please see [43].

Having run the circular flow simulation for $Re=4000$, we visualized its corresponding simulation data using VisIt [41], see Figure 24, which gives colormaps of vorticity, magnitude of velocity, velocity in the horizontal and vertical directions, and the finite-time Lyapunov exponent (FTLE). The FTLE can be used to find the rate of separation in the trajectories of two infinitesimally close parcels of fluid. Maxima in the FTLE have been used to determine approximations to the true Lagrangian Coherent Structures (LCSs). LCSs are used to determine distinct flow structures in the fluid [58,59,60,61] can be used to divide the fluid’s complex dynamics into distinct regions to better understand transport properties of flow [62,63,64]. True LCSs would require knowledge of the infinite-time Lyapunov exponent (ITLE); however, the FTLE can serve as a proxy to the ITLE. In this paper, we computed the forward-time FTLE field, whose maximal ridges give approximate LCSs corresponding to regions of repelling fluid trajectories and whose low values give rise to regions of attraction [61]. Thus, the FTLE can be used to help find regions of high fluid mixing [58,59,61,65,66].

Figure 24 also provides contours for each quantity. It also includes a depiction of the fluid velocity field, as scaled by the maximum in the entire domain’s magnitude of velocity, which also includes its associated streamlines. Recall that streamlines are curves that depict instantaneous tangent lines along direction of the flow velocity. They show the direction that a neutrally-buoyant particle would travel at a particular snapshot in time. These are different than contours (also known as level-sets or isolines), which are curves where a function has a constant value. This snapshot was taken at $t=24.0\mathrm{s}$, when the simulation ended.

As the velocity boundary conditions continually increased, the interior of the domain began to move clockwise, in the same direction as the flow at the boundaries (see Figure 1b). The colormaps and contours in each velocity-related panel reveal that flow speeds decrease towards the interior of the domain. Moreover, as all the boundary conditions are uniform, the vertical and horizontal velocity plots are identical under a 90 degree rotation. Also, the vorticity plot illustrates that vorticity is still low by the ending of the simulation in the center of the domain, even though flow across the whole domain is rotating clockwise. The FTLE plot suggests that the majority of fluid mixing occurs near the edge of the domain, rather than the interior, as higher values of FTLE suggest nearby fluid trajectories move away from each other at higher rates. This is due to larger spatial gradients in the velocity in these regions, where flow is moving at different speeds in different directions.

We also present data comparing simulations for $Re=400,1000$, and 4000 (for $\mu =2.5,1.0$, and $0.25\mathrm{kg}/(\mathrm{m}·\mathrm{s})$, respectively). Figure 25 illustrates colormaps of the magnitude of velocity (with corresponding contours) at various time-points during the simulation. Every colormap uses the same scaling in the figure. Higher velocities appear to creep into the interior of the domain quicker in the lower $Re$ cases. By $t=24\mathrm{s}$ it appears that the $Re=400$ and 1000 cases look almost identical, the $Re=4000$ case tells a different story - the interior of the domain is still moving considerable slower than the other two cases. This might seem counter-intuitive at first, as one might generally think that higher $Re$ tends to lead to more fluid motion, especially when we are lowering viscosity to increase $Re$. However, this is due to differences in time-scales for the diffusion of momentum through the fluid. The viscous diffusion time can be thought of as the time it takes for a fluid parcel to diffuse a particular distance on average. If $\tilde{t}$ is the diffusion time, $\nu =\mu /\rho$ is the kinematic viscosity, and $\tilde{l}$ is the mean-squared distance, then

$$\tilde{t}\sim \frac{\tilde{l}}{\nu }\sim \frac{1}{\nu },$$

if all of the simulation geometries are identical and only $\mu$ (and hence $\nu$) is varied. This concept is based off viewing diffusion as a random walk process [67]. Hence for the simulations presented here the time-scales vary between

$${\tilde{t}}_{Re=400}\sim \frac{1}{2.5/1000}=400\mathrm{s}\mathrm{and}{\tilde{t}}_{Re=4000}\sim \frac{1}{0.25/1000}=4000\mathrm{s}.$$

Therefore in the $Re=4000$ case, the diffusion time-scale is $10x$ longer than that of the $Re=400$ case, and thus the dynamics evolve much slower in the $Re=4000$ case! Moreover, this phenomena is shown in Figure 26, which illustrates the fluid velocity field along with its associated streamlines at the same snapshots in time.

This is more apparent if we compare flow profiles within the domain. Figure 27 compares the horizontal flow profile for $Re=400,1000$, and 4000 across the vertical mid-line of the domain. As mentioned above the dynamics in the $Re=4000$ case evolves about $4x$ and $10x$ slower than the $Re=1000$ and $Re=400$ cases, respectively. By $24.0\mathrm{s}$, the flow profiles in the $Re=400$ and $Re=1000$ cases look almost identical, while the flow profile in the $Re=4000$ resembles that of the the flow profile in the $Re=1000$ at ∼6.4 s—approximately a factor of 4 different in time. Furthermore Figure 28 shows the flow profiles for all three cases of $Re$ when the diffusion of momentum has reached the same depth. This occurs at different simulation times due do the variations in viscous diffusion time-scales.

Students may elect to try the following:

• Recreate the above results with the parameters listed in Table 4
• Change the computational geometry, e.g., rectangular instead of square
• Vary the input velocity along each side of the domain to make them non-uniform
• Vary the $Re$ by changing dynamic viscosity for suggestion (2) or (3)

4.3. Side-by-Side Voritices (via Spectral Method)

In this example we will use a FFT-based fluid solver in the software to run a simulation of multiple regions of vorticity interacting with one another. Note that the first example given here can be run by going into in the FFT_NS_Solver script selecting the ‘qtrs’ option (see line 47 in Figure 3b). The vorticity is initialized as in Figure 2a. We will also highlight the ‘half’ option in this section as well.

Recall that for this solver, an initial vorticity configuration is needed. For the simulations of four interacting vorticity regions, the $CW$ and $CCW$ vorticity were initialized as $\omega =-1$ and $1\mathrm{rad}/\mathrm{s}$, respectively. In other cases with only two interacting voricity regions, the strengths of each regions were varied between $-1,-0.5,0.5$ and $1\mathrm{rad}/\mathrm{s}$, as appropriate. For the simulations shown below with either 2 or 4 interacting vorticity regions, use the computational parameters listed in

Table 5. Numerical parameters for the 4 interacting vortices simulation.

Parameter	Variable	Units	Value
Domain Size (4 case)	$[L_x,L_y]$	m	$[1,1]$
Spatial Grid Resolution (4 case)	$[N_x,N_y]$		$[512,512]$
Spatial Grid Size	$dx=dy$	m	$L_x/N_x=L_y/N_y$
Domain Size (2 case)	$[L_x,L_y]$	m	$[1,0.5]$
Spatial Grid Resolution (4 case)	$[N_x,N_y]$		$[512,256]$
Spatial Grid Size	$dx=dy$	m	$L_x/N_x=L_y/N_y$
Time Step Size	$dt$	s	${10}^{-2}$
Total Simulation Time	T	s	5
Fluid Kinematic Viscosity	$\nu =\rho /\mu$	${\mathrm{m}}^2/\mathrm{s}$	$0.001$

.

Figure 29 provides the simulation data for the case with 4 interacting vorticity regions, as described above. It presents colormaps for vorticity, magnitude of velocity, horizontal velocity, vertical velocity, and the finite-time Lyapunov exponent (FTLE), which is used to determine approximations to the true Lagrangian Coherent Structures (LCS) to which can be used to distinguish fluid mixing regions [58,59,61,65,66]. Their corresponding contours are also given. Recall that contours are curves to which the value of a quantity is constant. Figure 29 also gives a snapshot of the velocity vector field. This data is from the last snapshot of the simulation at $t=5.0\mathrm{s}$.

The top two and bottom two vorticity regions are initialized the same, e.g., clockwise (CW) and counterclockwise (CCW), respectively, by this point in the simulation they are beginning to lean in the direction of rotation. Moreover, it can be seen that fluid is moving quickly horizontally through the top, middle, and bottom of the domain, as illustrated by the magnitude of velocity plot. When this is compared to the horizontal velocity plot, one observes that the fluid along the top and bottom of the domain both are moving left-to-right, while fluid moving horizontally through the middle of the domain is moving right-to-left. These directions align with the direction of the spinning vortices in each of these horizontally-aligned regions—top, middle, and bottom, which can be seen in the velocity vector field depiction.

The high values of FTLE illustrate regions of higher fluid mixing, as higher FTLE suggest faster rates of separation of close fluid blobs. In Figure 29 these are the areas between the top-two and bottom-two vorticity regions. Within these regions there is intense shearing among the vertical flow velocity, see the vertical velocity plot. Red regions in the vertical velocity plot indicate fluid moving upwards while blue corresponds to downward fluid motion. Moreover, from the nature of the FFT-solver, the periodic boundary conditions, illustrate that along vertical walls of the domain on opposite sides there are regions where fluid is moving the the opposite direction. Hence the large FTLE values on the left and ride side of the top-two and bottom-two vortices. There is significantly less mixing in the top, middle, and bottom horizontal strips across the domain because fluid is generally moving unidirectionally in these regions with lower spatial gradients.

One could elect to change the size or strength of a subset of these vorticity regions to simulate asymmetric vorticity interactions in this configuration. However, rather than study 4 interacting vorticity regions, we compared 2 interacting regions at a time, and varied the strength (magnitude of vorticity), size, and initial vorticity values (spinning-direction) between them, see Figure 30. Figure 30 provides comparison data of snapshots at $t=6.0\mathrm{s}$ in each simulation for different quantities—either vorticity with velocity vectors, magnitude of velocity with its corresponding contours, and FTLE with its corresponding contours. We use the language that a “strong” vorticity region has a $\left|\omega \right|=1\mathrm{rad}/\mathrm{s}$, while a “weaker” region has a vorticity of $\left|\omega \right|=0.5\mathrm{rad}/\mathrm{s}$. Each initial size of the vorticity regions had a radius of $0.15\mathrm{m}$ or $0.3\mathrm{m}$.

In general, one can notice the following:

• When same-spinning and size vorticity regions interact, more fluid mixing occurs directly between them (1st row)
• When oppositely-spinning and same size vortex regions interact, there is enhanced vertical fluid flow between the vorticity regions, but less mixing between them (2nd row)
• When oppositely-spinning vorticity regions of different strengths interact, the stronger region (whose magnitude of vorticity is greater), influences fluid flow to move in the same direction further away from it. There is still some enhanced vertical flow through the region between the regions (3rd row)
• If there are asymmetric region sizes and strengths, where the larger region is weaker (smaller magnitude of vorticity), the region with the larger vorticity magnitude may be able to keep the other region from influencing much of the flow around it, while imposing its flow direction on the weaker region (4th row)
• Two regions of vorticity with the same sign, but differing strengths, leads to enhanced fluid mixing between the vortices, while flow around the smaller region are significantly less than the other (5th row)
• If there are asymmetric region sizes but uniform strengths, enhanced vertical flow will be seen between the regions, although the larger vorticity region exerts more influence on flow patterns closer to the smaller region (6th row)
• If there are asymmetric region sizes and strengths, where the smaller region is weaker (smaller magnitude of vorticity), enhanced vertical flow will be seen between the regions, although the stronger vorticity region exerts more influence closer to the smaller region (7th row)

Students may elect to try the following:

• Recreate the above results for interacting side-by-side vorticity regions with the geometric and simulation parameters given in Table 5 or for different $\nu$ or computational domains.
• Add more vorticity regions or make the case with four asymmetric in terms of placement and/or size.
• Change the distance between the vorticity regions to observe how magnitude of velocity or FTLE contours change.

4.4. Evolution of Vorticity from an Initial Velocity Field

In this example using a spectral (FFT) solver, we begin with a snapshot of the velocity vector field from an independent CFD simulation, see Figure 31a,b. Note that this snapshot was taken from an open-source fluid-structure interaction example of a pulsing cartoon heart [36]. We used a stored time-point’s velocity field data, to which we extracted the horizontal and vertical velocity components, and then computed its associated vorticity at that particular time-step, see Figure 31c. We note that there are numerical errors in computing the vorticity, $\omega$, as a centered finite difference scheme [43] was used to compute each first-order partial derivative in calculating vorticity, i.e., $\omega =\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$. Moreover, since this fluid solver uses periodic boundary conditions, we computed the partial derivatives at each boundary using data from the opposite side of the computational domain. Note this example may be run by selecting the ‘jets’ option in the FFT_NS_Solver script (see line 47 in Figure 3b).

After having computed the vorticity from a velocity field, the simulation could begin. This simulation used the computational parameters listed in

Table 6. Numerical parameters for case of evolving vorticity from an initial velocity field.

Parameter	Variable	Units	Value
Domain Size	$[L_x,L_y]$	m	$[1,1]$
Spatial Grid Resolution	$[N_x,N_y]$		$[256,256]$
Spatial Grid Size	$dx=dy$	m	$L_x/N_x=L_y/N_y$
Time Step Size	$dt$	s	${10}^{-2}$
Total Simulation Time	T	s	50
Fluid Kinematic Viscosity	$\nu =\mu /\rho$	m${}^2$/s	$0.001$

. We further note that the velocity field had an original grid resolution of $[N_x,N_y]=[512,512]$, so we down-sampled it to a $[256,256]$ grid for the example shown here. Furthermore, as the original velocity field was computed from an independent fluid-structure simulation, we comment that there is no complex, moving boundary within the computational domain; however, at the initial point you can see the remnants of the coinciding heart structure, see Figure 31.

Figure 32 provides snapshots over the simulation to illustrate how the vorticity and magnitude of velocity evolve. Both Figure 32a,b illustrate the effect of periodic boundary conditions along each edge of the domain, e.g., vortices moving through the top or right side of the domain continue through the bottom or left side of the domain, respectively. Furthermore, in the snapshots provided, the overall vorticity and magnitude of velocity appears to dissipate within the domain. This is because there are no source terms (or explicit boundary conditions) providing additional flow (energy) into the system to induce more flow. Although the original velocity vector field that defined the initial vorticity of the system came from a simulation of pulsing hearts, without the heart moving and thus pushing on the fluid, the fluid will eventually stop moving and reach zero velocity due to the fluid’s viscous nature. Therefore initializing the vorticity out of a time-point from another simulation’s velocity field, is akin to studying the evolution of the fluid’s motion from a singular impulse in the flow.

Students may elect to try the following:

• Recreate the above results with the parameters listed in Table 6 or for different $\nu$.
• Take a velocity field from one of the Projection method examples, e.g., from Section 4.1 or Section 4.3, and use it to initialize the vorticity for a simulation using this spectral (FFT) method.
• Find how long it takes the simulation to reach a $1/10$ or $1/100$ of the original maximum magnitude of velocity.
• Discover how that time varies with changes in viscosity, $\nu$.

4.5. Flow Past One or More Cylinders (via Lattice Boltzmann)

Flow past objects is one of the most heavily studied problems in aerodynamics and thus is a standard example used in many fluid dynamics courses [68,69,70]. While many texts only consider the problem for inviscid flows and use potential flow methods to solve for exact solutions [69,70], here we include effects of viscosity to illustrate transitions to vortical flow, and upon doing so, a subset of the possible flow separation cases. Furthermore, flow past cylinders is also still an active area of contemporary research [71,72,73,74,75,76,77].

We performed simulations for both a single cylinder as well as multiple cylinders. The computational geometry for a single cylinder case is given in Figure 33. Figure 33a gives the boundary conditions considered for the problem for flow past a rigid cylinder contained within a rigid channel. The walls of the channel are rigid and so bounce-back boundary conditions are used. The ends of the channel use periodic boundary conditions, e.g., what goes out the right side, comes in through the left. Note that this is the first example discussed in which uses a mix of explicit boundary conditions and periodic boundary conditions. Moreover, the cylinder itself is modeled as a rigid structure and so bounce-back conditions are used on it. Figure 33b provides visual details how the inflow is modeled in the simulation. Every iteration, the horizontal velocity at the inflow is increased by an amount of $\Delta U_x$ to drive fluid flow towards the right.

All simulation parameters (fluid, grid, and geometry) for both cases involving either a single cylinder or multiple rigid cylinders are given in

Table 7. Numerical parameters for flow past one or more cylinders using the Lattice Boltzmann Method.

Parameter	Variable	Units	Value
Domain Size	$[L_x,L_y]$	m	$[2,0.5]$
Spatial Grid Resolution (1 Cylinder)	$[N_x,N_y]$		$[640,160]$
Spatial Grid Resolution (Multiple Cylinders)	$[N_x,N_y]$		$[512,128]$
Spatial Grid Size	$dx=dy$	m	$L_x/N_x=L_y/N_y$
‘Density’ Initialization	$\rho$		$0.01$
Relaxation Parameter	$\tau$		$0.53$
Total Simulation Steps (1 Cylinder)	N		56000
Total Simulation Steps (Multiple Cylinders)	N		5500
Incremental inflow velocity increase (1 cylinder)	$\Delta u_X$		$0.00125$
Incremental inflow velocity increase (Multiple Cylinders)	$\Delta U_x$		$0.01$
Radii (1 Cylinder)	r	$\%L_y$	$7.5\%L_y$
Radii (Multiple Cylinders)	r	$\%L_y$	$12.5\%L_y$

. Note that the cases with multiple cylinders simulate flow around cylinders with larger radii, use a larger value for $\Delta U_x$, and a different final number of steps and grid resolutions. Furthermore, due to these differences we are not comparing the case of one single cylinder to the case with multiple cylinders.

The simulation data obtained at step n = 49,600 is given in Figure 34. It presents colormaps (and corresponding contours) for vorticity, magnitude of velocity, horizontal velocity, vertical velocity, and the finite time Lyapunov exponent (FTLE), to which knowledge of the approximate Lagrangian Coherent Structures (LCS) can be extracted [58,59,61,65,66]. We attribute regions with high FTLE to regions of high fluid mixing, as higher FTLE suggests that nearby fluid parcels separate at a much faster rate than those in lower FTLE regions. From the vorticity panel, significant flow separation is seen as vortices are being shed off the cylinder and flow pushes past it left-to-right. Oppositely spinning vortices are shed off cylinder one after another. Moreover, we observe that there is significantly more horizontal fluid motion than vertical within the channel. Since flow is being pushed left-to-right from the inflow condition, it would require more energy for fluid to move vertically, than to simply…go with the flow. Figure 35 provides snapshots of vorticity (left) and FTLE (right) during the simulation to highlight how smooth flow developed into vortical flow patterns once the horizontal velocity increased past a threshold. Such threshold appears to have occurred around simulation step ∼47,000.

In particular, flow instabilities began to manifest by step 46,000, as seen by both the vorticity and FTLE snapshots. The majority of fluid mixing occurred in large contours behind the cylinder up to this point. Once flow separation occurs, some mixing occurs in the wake in small patches, but no geometrically long regions of higher fluid mixing are present, as seen by the contours. Note that in the last two panels, for steps 49,200 and 49,600, large FTLE-valued regions do not correspond to locations of shed vortices, rather higher FTLE valued regions appear between oppositely spinning vortices.

Finally we performed simulations with multiple cylinders in a channel. Each cylinder in these simulations had a radii 66% larger than the single cylinder shown above. We performed three separate simulations with different configurations of the cylinders, each based off a triangle formation. Figure 36 illustrates how flow evolved for each geometric configuration, in particular for (a) vorticity and (b) FTLE.

Vortices are observed shedding off each cylinder, possibly at an enhanced rate due to flow interactions with the other cylinders (see suggested activity below). In the symmetric configurations, e.g., the triangle-formation and vertical-line formation, flow symmetry is preserved, while in the case with only one trailing cylinder, flow asymmetries arise. FTLE plots at the last step shown (5000) illustrate different geometrical configurations not only can lead to different patterns of fluid mixing, but also enhanced mixing, even among cases with only two cylinders. The vertically-aligned cylinder case appears to elicit more downstream mixing in the wake than the case of asymmetrically placed cylinders; however, this may also be an artifact of overall geometry of the system, e.g., the size of the cylinders and the proximity of cylinders to the channel walls.

Students may elect to try the following:

• Recreate the above results with the parameters listed in Table 7
• Using the radius for the larger geometry case, run a simulation with only one cylinder present in the middle and compare vorticity and FTLE evolution plots.
• Vary $\tau$ (relaxation parameter relate to viscosity) to test the system’s sensitivity to $\tau$
• Change the number, placement, or size of each cylinder in the domain.
• Change the shape of the object in the channel, e.g., try an ellipse or airfoil-like shape

4.6. Flow Past a Porous Cylinder (via Lattice Boltzmann)

While flow around solid objects is a popular problem in aerodynamics, flow around porous objects has only recently begun to be investigated within in the past few decades, using either high fidelity numerical simulation or experiments [78,79,80,81,82]. Here, we present an example of flow past a porous cylinder using the LBM. The cylinder geometry, computational setup, and inflow/boundary conditions are identical to the single cylinder case of Section 4.5 (see Table 7). The only difference being that the cylinder geometry itself is porous.

The porosity, or void fraction ($VF$), of the cylinder was chosen to be $VF\in \{0\%,10\%,25\%,50\%,62.5\%\}$. To model the void fraction, each grid cell inside the cylinder was given a random number between $[0,1]$ and then depending on the specific $VF$, if that grid cell’s randomly assigned number was above the $VF$ in decimal form, it was assigned as a solid boundary. For example, for $VF=25\%$, each grid cell with a randomly assigned number greater than $0.25$ was declared a solid boundary point. Figure 37 depicts the porous cylinders considered for each $VF$ case studied below.

Although, the cylinder is no longer uniformly solid, in this example we will not focus our efforts to study flow through the cylinder itself. For example, in Figure 37’s $VF=10\%$ case, there are holes in the interior of the cylinder domain that are not connected to the fluid region outside of the cylinder. Those void regions will not influence the flow outside of the cylinder, as they are blocked from outside fluid penetrating them. However, if you wanted to study the fluid dynamics through a porous structure, you could design a specific void network structure within the cylinder or in another desired geometry. What we will emphasize here is that the porous regions connected to the outside of the cylinder cause asymmetries to develop in the overall flow pattern at an accelerated rate as compared to the non-porous (solid) case. Such asymmetries then lead to quicker vortex shedding.

Figure 38 shows the evolution of vortex shedding in cases of differing porosities (void fractions) of $\{0\%,10\%,25\%,50\%,62.5\%\}$. Every case shows that more porosity leads to faster development of vortex shedding. Furthermore, even the case with $VF=10\%$ has pronounced vortices being shed before the $VF=0\%$ case even shows any signs of significant flow asymmetry. The increased porosity accelerates flow asymmetries to manifest, leading to vortex shedding occurring earlier on. This example highlights how small perturbations in fluid dynamics problems (here small differences in geometric structure) can lead to the system having significantly different time-scales for flow structures to develop or for a bifurcation in the resulting dynamics to emerge.

Students may elect to try the following:

• Recreate the above results with the parameters listed in Table 7
• Create the accompanying FTLE plots and discuss how fluid mixing changes due to the porous cylinder
• Vary $\tau$ (relaxation parameter relate to viscosity) to test the system’s sensitivity to $\tau$
• Increase the number of porous cylinder in the domain and vary their placement or size
• Change the shape of the porous object in the channel, e.g., try an ellipse or airfoil-like shape
• Create a specific porous network structure through the cylinder and test how flows directly through a porous cylinder may affect vortex shedding

5. Discussion

In this work we provide software that was developed for the specific purpose to make CFD accessible to undergraduate (and graduate) students in order to provide them an opportunity to perform traditional and contemporary CFD simulations as a valuable learning experience. To that end all the fluid solvers contained within were written in both MATLAB and Python, two languages that are commonly familiar to most science and engineering students [24,25,26,27,28,31]. Students also have the opportunity to modify existing examples or create their own examples within the software to test hypothesis or for stimulate further scientific curiosity. To that extent, we provided an overview of how the code is structured (Section 3) and offered interesting variations to each of the examples that we showcased in this paper (Section 4).

Not only can students run CFD simulations, they also have the chance the practice the art of scientific visualization, and use open-source software that is used by many researchers (VisIt [41] or Paraview [42]), to visualize the corresponding data. These software packages are both open-source and have easy-to-use and learn graphical user interfaces (GUIs) that streamline the visualization process. They can also be used to perform data analysis as well, which further reduces the computational learning curve for students to be able probe the data produced. We provided step-by-step guides on how to use VisIt for these tasks in Section 3.

By performing their own CFD experiments, fluid dynamics students have the opportunity to vary parameters (e.g., fluid scale ($Re$), boundary conditions, geometry, etc.) of a given system and observe the resulting dynamics. More specifically, they are able to decrypt how the dynamics may vary across parameter regimes, both qualitatively and quantitatively. This grants them the ability to challenge their own developing intuition of the resulting flow physics to further accelerate their learning. They also gain more computer experience and witness first-hand some of the advantages that complex computer simulations offer.

If students wish to go a step further and study the underlying mathematical structure of the code, we provide insightful comments everywhere possible within each fluid solver script. To complement this, we provided a detailed mathematical overviews of each fluid solver used here in Appendix B. While it is not our mission here to teach students how implement their own numerical fluid solvers, we believe this work contributes to providing a foundation for capturing the importance (and utility) of different fluid solver schemes. For each solver that was introduced for particular applications in Section 2.1, Section 2.2 and Section 2.3, we give a high-level overview of the method. If students wish to continue learning about numerical methods for solving the fluid equations or numerical methods for partial differential equations in greater depth, we suggest Barba and G. Forsyth’s CFD Python: the 12 steps to Navier-Stokes equations [33], L.A. Barba and O. Mesnard’s Aero Python: classical aerodynamics of potential flow using Python [34], both of which use Jupyter Notebooks [38], or Pawar and San’s [39] modules for developing fluid solvers in the Julia programming language [40].

In this day and age, as the importance of computer literacy increases rapidly [83,84,85,86], integrating more hands-on computer-based activities into course structures will be of critical importance. However, students also seeing successful execution of code is paramount [85]. Here we provide gateway software for students to dive into the world of CFD, in familiar programming environments. While proprietary software offers immense advantages to running simulations and analyzing data, students may see a disconnect between the programming knowledge they’ve acquired and such polished software. We hope to provide them a unique opportunity to see fluid solvers written in familiar, non-intimidating environments that they may be able to tweak and modify comfortably. Our hope is to inspire further ownership of their learning and to stimulate more interest in (lucrative, transferable) computer programming skills.