This article describes how Mathematica can be used to develop an understanding of the basic steps involved in solving Navier–Stokes equations using a finite-volume approach for incompressible steady-state flow. The main aim is to let students follow from a mathematical description of a given problem through to the method of solution in a transparent way. The well-known “driven cavity” problem is used as the problem for testing the coding, and the Navier–Stokes equations are solved in vorticity-streamfunction form. Building on what the students were familiar with from a previous course, the solution algorithm for the vorticity-streamfunction equations chosen was a relaxation procedure. However, this approach converges very slowly, so another method using matrix and linear algebra concepts was also introduced to emphasize the need for efficient and optimized code.

Introduction

Mathematica is used to help with an initial understanding of the process of solving Navier–Stokes equations. For 2D incompressible flows, it is possible to recast the Navier–Stokes equations in an alternative form in terms of the streamfunction and the vorticity. In many applications, the vorticity-streamfunction form of the Navier–Stokes equations provides better insight into the physical mechanisms driving the flow than the “primitive variable” formulation in terms of the mean velocities , , and pressure . The streamfunction and vorticity formulation is also useful for numerical work, since it avoids some problems resulting from the discretization of the continuity equation.

The streamfunction is defined as

(1)

where the integral has to be evaluated along a curve from the arbitrary but fixed point to point , is the velocity vector, and is the unit normal on the curve from to ; see Figure 1. We regard as a function of the location of point .

Figure 1. Sketch illustrating the definitions of a streamfunction.

Figure 1 shows that is equal to the component of the velocity that crosses . Therefore represents the volume flux (per unit depth in the direction) through . Evaluating along two different paths and invoking the integral form of the incompressibility constraint shows that is path independent; that is, its value only depends on the locations of the points and . Changing the position of point only changes by a constant. It turns out that for all applications such changes are irrelevant. It is therefore common to suppress the explicit reference to . Hence, we regard as a function of the spatial coordinates only; that is, . Streamlines are lines that are everywhere tangential to the velocity field, that is, , where is the unit normal to the streamline. Hence the streamfunction is constant along streamlines. Note that stationary impermeable boundaries are also characterized by , where is the unit normal on the boundary. Therefore, is also constant along such boundaries. Invoking the integral incompressibility constraint for an infinitesimally small triangle shows that is related to the two Cartesian velocity components and via

(2)

Flows that are specified by a streamfunction automatically satisfy the continuity equation, since

(3)

For 2D flows, the vorticity vector only has one nonzero component (in the direction); that is, , where

(4)

Using the definition of the velocities in terms of the streamfunction shows that

(5)

where is the 2D Laplace operator.

The understanding of the steps involved in solving the Navier–Stokes equations using the vorticity-streamfunction form is one of the topics used in a third-year undergraduate course on computational fluid mechanics, solely for students majoring in mechanical engineering. The use of Mathematica makes the assumption that the students are familiar with the package, as it generally takes a good deal of exposure to Mathematica to become comfortable using it at the level required here [1]. The students taking the computational fluid mechanics course are indeed very familiar with Mathematica, as the computer algebra system is used during year one in the modules Calculus and Applications and Vector Calculus, and during year two in the module Numerical Methods for Engineering. In addition, there are many notes, explanations, and examples on the in-house Moodle open-source learning platform. Similar work to this can be found in Fearn [2], while background reading on computational fluid dynamics and vorticity may be found in Ferziger and Perić [3] and Chorin [4], respectively.

The aim here is to give a continuous and comprehensive process involving the mathematics of one formulation of the Navier–Stokes equations and a solution using Mathematica. In this way, the students can actually see the development of the mathematical description linked to a programming environment solution process so often hidden in commercial code used for training CFD students.

Navier–Stokes Equations in Vorticity-Streamfunction Form

The Navier–Stokes equations for incompressible steady-state flow in vorticity-streamfunction form are

Advection-Diffusion Equation

(6)

Elliptic Equation

(7)

where is the Reynolds number. For a given problem, boundary conditions need to be specified to solve equations (6) and (7). The problem chosen here introduces the student to a classical computational fluid dynamics (CFD) case, namely the lid-driven cavity flow [5]. This flow is commonly used to test, for example, a novel method of discretization of the equations or new computer programming, as the resulting flow is well known from experiments. Consider a rectangular box as shown in Figure 2, where a lid is allowed to move in the horizontal plane from left to right. When the lid is not moving, the fluid in the box is stationary, whereas when the lid is moving, the fluid circulates inside the box. Here the boundary conditions needed for solution are summarized in terms of vorticity and streamfunction. As all four walls of the cavity touch, the streamfuction must be equal for all four walls, as indicated by the in Figure 2. The streamfuction is constant at the walls, as its gradient is velocity, which is zero relative to a given wall (no-slip condition), as indicated in the outer boxes also shown in Figure 2. The vorticity boundary conditions for each wall are also shown in each of these outer boxes and were derived from the streamfunction.

Figure 2. Summary of boundary conditions in terms of vorticity and streamfunction for cavity with moving lid.

Equations in Dimensionless Form

It is convenient from a numerical point of view to make equations (6) and (7) nondimensional. If and are the height and width of the cavity, respectively, and is a reference velocity, then

.

Equations (6) and (7) become

(8)

(9)

defined on the domain . The boundary equations are defined by

(10)

Discretization and Solution Algorithm

Discretization of Transport Equations

Here is an iterative solution for solving the incompressible steady 2D Navier–Stokes equations. The method is based on Burggraf’s proposal [6]. So as to include the boundary nodes of the cavity, central differencing is used to discretize equations (8) and (9). The number of nodes in the and directions is set at and , respectively. The mesh is a regular Cartesian grid with the nodes equally spaced at a distance , as illustrated in Figure 3.

Figure 3. Typical Cartesian mesh used for the lid-driven cavity flow.

To discretize equations (8) and (9), the following general finite central-difference approximations were introduced for spatial dimensions, obtained by adding or subtracting one Taylor series from another and invoking the intermediate value theorem. As can be seen from equation (11), the order of accuracy for both the first and second derivatives is .

(11)

Using the approximations from equation (11), the discretized advection-diffusion equation (equation (8)) at the node is

(12)

The elliptic equation (equation (9)) at the node is

(13)

The equations are applied to the internal nodes of the Cartesian mesh; that is, and .

Solution Algorithm

The solution method uses residual functions; that is, if the values of and are exact on the nodes spanned by the residual functions and , then

(14)

where

(15)

(16)

The following fixed-point iterative procedure, based on the Gauss–Seidel scheme, is then constructed:

(17)

where is a relaxation parameter lying in the range , and and refer to the respective iterations. In this article, the actual value of depends on and can be obtained by numerical experimentation. The use of this relaxation parameter is to improve the stability of a computation, particularly in solving steady-state problems. It works by limiting, when necessary, the amount that a variable changes from one iteration to the next. The optimum choice of is one that is small enough to ensure stable computation but large enough to move the iterative process forward quickly.

The boundary conditions now need to be determined. For the streamfunction from equation (11),

(18)

For vorticity, the following needs to be considered for, say, the left wall of the cavity with one node outside the computational domain:

(19)

The value of can be accounted for using the no-slip boundary condition on the cavity wall; that is,

(20)

Using central differences, this condition can be written as

(21)

where again the order of accuracy is .

This finding, together with , gives

(22)

Similarly, for the other walls of the cavity the boundary conditions for vorticity can be deduced from

(23)

to give

(24)

where and are defined in Figure 3. Note that the final boundary condition shown in equation (24) is for the moving lid, and the extra term is due to the tangential velocity being nonzero.

Use of Mathematica

This section starts with the definition of mesh size and grid spacing; for this problem the mesh spacing is set equal in the and directions. The parameters Nx and Ny denote the mesh size. Initially, ω and ψ are set to zero for all nodes except at the top wall, where vorticity is not zero. The Reynolds number is set at 100, the relaxation parameter p is set to 1, and the aspect ratio γ is set at 1. A maximum allowable residual value e is set, and nested loops are used to execute the iterative algorithm. At every iteration step k, if the maximum value of the absolute value of is less than e, the calculations are halted. The iterative loops are wrapped with the function Timing to give an estimate of the time taken to do the calculations.

Geometry, Mesh Parameters, Initial and Boundary Conditions

Transport Equations

Iterative Algorithm

Streamfunction Results

Vorticity Results

Extraction of the centerline velocities is also instructive. First of all, students are given more experience viewing velocity profiles, and second, there are ample calculations and measurements in the literature [7, 8, 9] for comparison with the results obtained here. The centerline velocities (in the direction) and (in the direction) were derived from the streamfunction values using the equations

(25)

(26)

The velocity distributions for and in the following figure are drawn along the vertical and horizontal centerlines, respectively.

Importantly, the students were made aware that the profiles they calculated must be compared, preferably with experimental measurements, to test the validity of the calculation technique. As can be seen from the preceding figure, the calculations for the horizontal centerline velocities compare very favorably with those reported in the literature [7, 9] for .

A More Efficient and Optimized Solution Method

It can be clearly seen that the relaxation procedure, though reasonably easy for students to follow and therefore educationally productive, is indeed very slow to converge. However, students must be aware that code should be written efficiently and optimized, making use of all available concepts and procedures. Therefore, as part of their understanding of solving Navier–Stokes equations in the best way, they were also introduced to the following Mathematica code. The preceding exercise was a small-scale problem using a fairly primitive method of solution, and even for this simple problem, the convergence time is prohibitive. Therefore, for large numerical computations and those with more complex geometry, it is important to use a solution method that will run more efficiently. This can be achieved using linear algebra and matrices and highly optimized functionalities built into Mathematica, namely NDSolve and LinearSolve. The following program solves the same problem as that described in Figure 2, except the Reynolds number is set at 400 and the computational grid is slightly more dense. The time to solve this same problem was reduced by over an order of magnitude. It should be understood that what the code is solving is the unsteady 2D Navier–Stokes equations. Each iteration is a time step from impulse-starting conditions to steady state, and for the lid-driven cavity flow, the unsteady solution converges to the steady state. The animation shown is for the streamfunction with a Reynolds number of 400.

Optimized Coding

Having established the method of solution, the student would then experiment by varying the grid size and seeking a solution independent of the grid, which is very important in CFD calculations. Results would also be obtained for different aspect ratios of the container. Again, where possible, calculations must be compared with experimental measurements to establish validity. As can be seen, the calculated horizontal velocity along the vertical centerline for was in good agreement with experimental results [7, 9].

Summary

This article outlines a well-defined sequence of steps needed to solve the Navier–Stokes equations cast in the vorticity-streamfunction form. The sequence is integrated with the use of Mathematica at appropriate stages to take care of the tedious computations, and hence to allow the students to concentrate on the overall details of the solution process. In addition, it is now important to integrate computer technology so as to complete lectures and theory. This has the advantages of helping with the computations, aiding presentation for reports and analysis, and motivating students. Incorporating Mathematica also takes away the “black-box” approach so often being used by students with full CFD commercial codes, which give no real understanding of the numerics involved. The idea of efficient and optimized coding was also introduced.

References

[1]	S. Pomeranz, “Using a Computer Algebra System to Teach the Finite Element Method,” International Journal of Engineering Education, 16(4), 2000 pp. 362–368. www.ijee.ie/articles/Vol16-4/IJEE1162.pdf.
[2]	R. L. Fearn, “Airfoil Aerodynamics Using Panel Methods,” The Mathematica Journal, 10(4), 2008 pp. 725–739. www.mathematica-journal.com/2008/11/airfoil-aerodynamics-using-panel-methods.
[3]	J. H. Ferziger and M. Perić, Computational Methods for Fluid Dynamics, 3rd ed., Berlin: Springer, 2002.
[4]	A. J. Chorin, Vorticity and Turbulence (Applied Mathematical Sciences, Vol. 103), New York: Springer-Verlag, 1994.
[5]	J. D. Bozeman and C. Dalton, “Numerical Study of Viscous Flow in a Cavity,” Journal of Computational Physics, 12(3), 1973 pp. 348–363. doi:10.1016/0021-9991(73)90157-5.
[6]	O. R. Burggraf, “Analytical and Numerical Studies of the Structure of Steady Separated Flows,” Journal of Fluid Mechanics, 24(1), 1966 pp. 113–151. doi:10.1017/S0022112066000545.
[7]	U. Gia, K. N. Ghia, and C. T. Shin, “High-Re Solutions for Incompressible Flow Using the Navier–Stokes Equations and a Multigrid Method,” Journal of Computational Physics, 48(3), 1982 pp. 387–411. doi:10.1016/0021-9991(82)90058-4.
[8]	W. F. Spotz, “Accuracy and Performance of Numerical Wall Boundary Conditions for Steady, 2D, Incompressible Streamfunction Vorticity,” International Journal for Numerical Methods in Fluids, 28(4), 1998 pp. 737–757. onlinelibrary.wiley.com/doi/10.1002/%28SICI %291097-0363%2819980930 %2928:4%3 C737::AID-FLD744 %3 E3 .0.CO;2-L/abstract.
[9]	D. C. Wan, Y. C. Zhou, and G. W. Wei, “Numerical Solution of Incompressible Flows by Discrete Singular Convolution,” International Journal for Numerical Methods in Fluids, 38(8), 2002 pp. 789–810. doi:10.1002/fld.253.
D. Adair and M. Jaeger, “Developing an Understanding of the Steps Involved in Solving Navier–Stokes Equations,” The Mathematica Journal, 2015. dx.doi.org/doi:10.3888/tmj.17-8.

About the Authors

Desmond Adair is a professor of mechanical engineering in the School of Engineering, Nazarbayev University, Astana, Republic of Kazakhstan. His recent research interests include investigations of airborne pollution for both passive and reacting flows, and developing engineering mathematics by the incorporation of computer algebra systems.

Martin Jaeger is an associate professor of civil engineering and head of department in the School of Engineering, Australian College of Kuwait, Misref, Kuwait. His recent research interests include construction management and total quality, as well as developing strategies for engineering education.

Desmond Adair
School of Engineering
Nazarbayev University
53 Kabanbay batyr Ave.
Astana, 010000, Republic of Kazakhstan
dadair@nu.edu.kz

Martin Jaeger
School of Engineering
Australian College of Kuwait
Al Aqsa Mosque Street
Misref, Kuwait City, Kuwait
m.jaeger@ack.edu.kw