## Using *Mathematica* to Simulate and Visualize Fluid Flow in a Box

July 9, 2013 — Paritosh Mokhasi, Mathematica Algorithm R&D

The motion of fluid flow has captured the interest of philosophers and scientists for a long time. Leonardo da Vinci made several sketches of the motion of fluid and made a number of observations about how water and air behave. He often observed that water had a swirling motion, sometimes big and sometimes small, as shown in the sketch below.

We would now call such swirling motions vortices, and we have a systematic way of understanding the behavior of fluids through the Navier–Stokes equations. Let’s first start with understanding these equations.

We have all seen special effects and visualizations of fluids in movies. The special effects are meant to give the viewer a realistic look at the motion of water or air. Since the focus in movies is on aesthetics, the flow may not be accurate. So, how does one go about writing code that will simulate fluid flow? In general, the answer is quite complicated, but with a few assumptions can be made more tractable.

Assuming that the simulation will take place under normal conditions, that is, at a reasonable temperature and at sea level, the equations that govern the motion of a fluid are:

These represent a time-dependent system of partial differential equations (or PDEs, for short) that can be solved with appropriate boundary conditions. *u* (m/s) , *v* (m/s), and *w* (m/s) represent the different components of the velocity; *p* (N/m^{2}) represents the pressure; *ν* (m^{2}/s) is the viscosity of the fluid; and *ρ* (kg/m^{3}) is the density of the fluid. Fluids obeying this system of PDEs are called incompressible fluids.

Solving the above system in general is difficult. Many people might decide at this point that the problem is much too complicated, and they don’t have the time or expertise to write specialized code to do this. This may have been true in the past, but with *Mathematica*, the simulation can be done using a number of built-in functions and visualization in a couple of lines.

In this blog post, we”ll look at some techniques to help in simulating fluid flow and we’ll visualize the motion of the flow.

Even though the problem seems intimidating, a lot of information can be gathered by just looking at the system. Let’s use an understandable example. Consider a square box where the top lid is allowed to move in the horizontal plane. When the top lid is not moving, the air (or fluid) inside is stationary. However, when the lid starts moving, the motion of the lid makes the fluid circulate inside the box. Now, even without understanding how the flow will move, we can make a qualitative guess about what might be happening inside the box. Here is a sketch of that.

The top lid drags the fluid below it, which generates a swirling motion. This initial swirling motion can instigate its own swirling motion, like a big gear turning smaller gears.

For the flow in a box, the Navier–Stokes equations can be simplified. First, we assume that the flow is not changing in time, and second, that it is not changing in the *z* direction. This simplifies the equations to:

This leads to three unknowns and three equations. If the velocities, pressure, and coordinates are non-dimensionalized as

where a characteristic velocity *U*_{0}(m/s) and a characteristic length *L* (m) are used, the equations become

where Re is called the Reynolds number and is defined as

As it turns out, the behavior of a flow can be described, to a large extent, by just this one number. The non-dimensionalization helps generalize the problem. So now the box can be made as large or as small as desired, and the top lid can be made to move fast or slow simply by changing the Reynolds number. This number can also be used to describe whether the fluid inside the box is water, air, or any other substance of your choice.

The Reynolds number is a very useful quantity in the field of fluid mechanics. This number can also be used to compare and characterize different flows. What does this mean? Suppose an engineer wanted to test the performance of a Boeing 747. It would be highly impractical (not to mention expensive) to run tests on the real airplane. The engineer can instead perform the tests on a model airplane with the Reynolds number of the model and the actual plane kept the same. This is the basis on which all wind-tunnel experiments are performed.

Now, back to the problem at hand. How can one solve the above system? As of now, there are no known analytical solutions, so the only choice is to obtain a solution using numerical methods.

The challenge in obtaining a numerical solution is to figure out the treatment of the derivatives. The general idea/algorithm is as follows:

1) Take the domain and discretize it. This means that the solution *u*(*x*) is computed at discrete points *x _{i}*, and the solution is to be found at

*u*(

*x*).

_{i}2) The derivatives are expressed as some linear combination of *u*(*x _{i}*).

3) The final system is changed from a system of differential equations to a system of algebraic equations for *u*(*x _{i}*). These can be solved using one of many root finding algorithms.

Let’s follow the algorithm outlined above. The first step is to discretize the domain. This means that we have to define a set of grid points on which the solution will be found. This is done as:

At each of these grid points, the velocities *u*(*x _{i}*,

*y*) and

_{i}*v*(

*x*,

_{i}*y*) and pressure

_{i}*p*(

*x*,

_{i}*y*) have to be computed. This means that each point contains three unknowns. We can give each unknown a name.

_{i}The second step in the algorithm is to find an approximation to the derivatives. The derivatives need to be approximated as:

This form of approximation is called finite difference approximation. The coefficients can be computed using the traditional Taylor series. `NDSolve` has a utility function called `NDSolve`FiniteDifferenceDerivative` that automatically computes these coefficients. Here is an example:

The function `NDSolve`FiniteDifferenceDerivative` can also be used to get the result in a matrix form. Having a matrix representation would make it easier to code. As an example, let’s get the first-order derivative for Sin(*x*) where *x* ∈(0,1).

So, with the help of `NDSolve`FiniteDifferenceDerivative`, the system of PDEs describing the flow can be converted to a system of nonlinear algebraic equations. For discretization, we need the differentiation matrices for the first and second derivatives in *x* and *y*. The accuracy of the finite difference approximation can be specified by the option “`DifferenceOrder`“. A value of 4 would mean that the approximation is 4th-order accurate.

With the variables defined and the matrices computed, the nonlinear algebraic equations are constructed as:

Notice that the Reynolds number is taken to be 100 (this can, of course, be changed).

There is one additional step (which turns out to be an important step) to account for the boundary conditions. To do this, we need to find the grid points that lie on the boundary AND the equation numbers associated with them. This can done as:

Now that the boundary positions are known, the following boundary conditions need to be applied:

Since the equation numbers associated with the boundaries are known, incorporating the boundary conditions involves a replacement of the discretized governing equations with the boundary equations.

With the discretization complete, we can proceed to solving the system. To do this, all the discretized equations should be joined into one large system of equations. This means that the final system to be solved will be 3*n* where *n* is the number of grid points. The 3 comes from the fact that we are solving three unknown variables at each grid point. So for a 100 x 100 grid, we will simultaneously solve 30,000 equations.

The function `FindRoot` can be used to solve this system.

Remember that the velocities are computed at discrete points. Using the function `Interpolation`, the discrete solution can be converted into a black-box continuous function.

We can visualize the interpolation using `LineIntegralConvolutionPlot`, which generates very attractive flow visualizations.

>

We can construct a number of helper/utility functions in order to put the code in a reusable form.

With the help of the utility functions, the main solver function is relatively compact.

The main solver function now consists of five parts. The functions of each of these parts are described as comments. Note that the main solver function “`DrivenCavitySolver`” takes in four arguments. These arguments are: Reynolds number, the domain, grid points in each direction, and the differentiation order. Having it in this form makes it easy to play with code. For example, the aspect ratio of the box can be changed to make it narrow or wide, and we can see what happens to the flow as a result of that change.

With the function now in hand, a number of interesting questions can be answered and different tests can be performed. The most obvious question that one can ask is: “What happens to the flow when the Reynolds number is increased?” Here are some results that easily answer this question through simulation and visualization.

Comparing the four flows, the common characteristic is a large central vortex dominating the flow. The differences between the four flows are seen in the secondary vortices. As the Reynolds number increases, the bottom-left vortex starts increasing in size, and eventually, when the Reynolds number is high enough, a third vortex starts forming at the top-left of the box (look at the fourth plot). We can therefore suspect the presence of more vortices as the Reynolds number is increased.

What we are viewing is a steady-state flow, which means that the flow does not change from one point in time to the next. This doesn’t mean that the flow is not moving! To visualize the actual motion of the flow, the box can be seeded with mass-less particles, and the particles are tracked in time as they are carried by the flow. This again can be done very easily using `NDSolve`.

The first step is to generate a flow. Let’s choose the flow at Reynolds number 2,000. This simulation will take a couple of minutes.

The box is seeded with 400 particles equally spaced throughout the box. The results of tracking the trajectories of each of these particles are called pathlines.

Using `NDSolve`, each individual particle can be tracked as it is carried by the flow. Let’s select a time window of 10 time-units to see how the particles move during that time period.

The evolution of the particles can now be visualized using the function `ListAnimate`.

Another idea for visualization is to release a continuous stream of dye at one location in space for a short period of time. The location can be determined by the experimenter. We will choose the location to be near the bottom-right of the box. A simple modification to the equations can achieve this effect.

Let’s experiment a little bit more. What happens when the aspect ratio is changed? If the aspect ratio is 2, then the movement of the primary vortex generates a second major vortex below it. As the strength of the primary vortex increases, it will tend to push the secondary vortex toward the left wall. Notice that if there is enough strength in this secondary vortex, it will generate its own vortex—the vortex at the bottom-left of the box.

What if the box is now made wide instead of narrow?

If you were to google “2D Navier Stokes Solver,” you would find hundreds if not thousands of papers and many algorithms. This can be overwhelming. My motivation for writing this short code was to see what results (if any) I could get if I were to take the raw equations and just solve them. With all the functions that *Mathematica* has to offer, and with a little numerical knowledge, I was able to get a prototype working pretty quickly. More importantly, *Mathematica* allowed me to explore the dynamics of the flow rather than sitting and coding for hours. I will be the first to admit that it is not optimal code, and many more features can be added to it. However, the ability to get something up and running quickly was empowering.

Download this post as a Computable Document Format (CDF) file.

## 41 Comments

I never, ever comment on blogs but let me say this post was just awesome. Wolfram Research needs to push Mathematica as the language of choice for numerical analysis. Many researchers in academia and industry view Mathematica as purely a symbolic analysis tool.

Right now there needs to be a strong marketing push to people that write these types of codes in MATLAB/C/what have you and show how Mathematica can be used to write CFD/FEA codes or to post-process what other analysis tools generate.

Thanks for the comments. I do plan on doing more posts highlighting the numerical side of Mathematica.

This is a great post. I agree that Mathematica needs to market its Numerical capabilities. This means focusing on the Compiler and highlighting the compilation to C capabilities. It would also be great if more mathematica commands can be compiled down to C. For example, Random numbers from different distributions and the Linear Algebra stuff (Linear Solve, Inverse, CholeskyDecomposition) etc could be made compilable. At the minimum, the functionality that is available from the Rmath standalone library of R could be compilable to expand Mathematica’s statistical capabilities. CUDALink could also be emphasized.

Please post on 2d transient heat conduction

You can solve transient heat conduction directly using NDSolve. Here is an example:

sol = NDSolveValue[{D[u[x, y, t], t] == 0.1*(D[u[x, y, t], x, x] + D[u[x, y, t], y, y]),

u[x, y, 0] == Exp[-8*(x^2 + y^2)], u[x, -1, t] == u[x, 1, t] == u[-1, y, t] == u[1, y, t] == 0},

u, {x, -1, 1}, {y, -1, 1}, {t, 0, 1}];

Animate[Plot3D @@ List[sol[x, y, t], {x, -1, 1}, {y, -1, 1}, PlotRange -> {{-1, 1}, {-1, 1}, {0, 1}},

PerformanceGoal -> “Quality”], {t, 0, 1}]

Thank you Sir

Hello, in the third of the Navier-Stokes there is a typo. You have -PartialD(p,x) when it should be -PartialD(p,z). Otherwise, very great and useful blog post!

Thanks for catching the typo. I corrected the equation.

Amazing!

Very interesting post!

Can’t wait for the follow up: transient (time dependent), 3D, inlet/outlet boundary, free surface (VOF), etc.

And also the other approaches: SPH (CPU vs GPU), Boltzmann, etc.

I also never post, but this one is amazing. Thanks!!!!

This is great piece of work! I would love to see more posts like this one. Thanks!

Ah – you are absolutely right. Thanks for catching that. It has been changed now.

Great blog post. Please do more.

I know that you couldn’t really get into analysis of the necessary grid size and time steps, but it would have been nice to indicate this.

excellet work i am love in it i have allready done the same story on csharp and was looking for it in mathematica it is awesome and simple as compared to csharp

.Great Blog Post………Very good tutorial….It is very exciting to solve a system of partial differential equations ..and specially NS equations (mother of fluid dynamics) !!!!….

Please post on….natural convection in a rectangular cavity..2D

I’m very interesting about your work.

I need your help for solving a doublet + uniform flow using the Navier-Stokes equations in 2D.

How will you impose the boundary conditions and simplify the equations ?

Thanks in advance.

Doublet+Uniform flow would be a potential flow. Potential flows have closed form solutions which can directly be used. There is no need to perform a numerical simulation. If on the other hand you are looking to solve viscous flow over a 2D cylinder then it will be a lot more complicated due to the the irregular geometry.

Isn’t the difficulty only how to impose the boundary conditions .. ?

and how can I solve the problem of the irregular geometry..

If you are looking to simulate viscous incompressible flow (i.e., solving Navier-Stokes eqautions) in the presence of irregular geometry, then we have to use an entirely different set of numerical tools to accommodate the domain AND more importantly the boundary conditions. In such cases, almost always, we end up using finite element methods to solve the problem. However, in my opinion, using immersed boundary methods is also equally competitive. Using such methods, boundary conditions just become part of the system of equations. I hope to put another blog together that discusses these immersed boundary methods in the very near future.

If, however, you are looking to solve inviscid potential flows, then using vortex panel/source methods are the most popular ones. These can handle any arbitrary geometry and you can enforce just about any boundary condition you like. The theory for such things is very well established and I think just about any book on Potential Flows would give you an algorithm for it.

Hope this helps.

How might this procedure be modified to estimate turbulence over an airfoil? My goal is to provide a resonable estimate of turbulence (e.g. Reynolds values) of airflow over a rotorcraft blade smooth leading edge and trailing surface as opposed one that surface anomolies at various distances from the center line, sizes and surfaces as the result of impact damage and other protective/sarificial surfaces.

Any help to a Mathematica novice, and one who has forgotten 90% of his fluid dynamics would be greatly appreciated.

Simulating flow over an airfoil would not be possible with the current code. It is designed specifically for flow in rectangular domains. Further, if the objective is to study turbulence, then the 2D NS equations would not help because turbulence cannot be generated by 2D flow (no momentum source term is available to generate vorticity). You would have to rely on some modelling methods like RANS or some variation of it to do your simulations.

Interesting post. What are the boundary conditions for P ?

Ah – I tried to avoid the discussion on pressure boundary conditions :). The pressure boundary conditions become apparent only after generating the pressure-Poisson equation (PPE), at which point, for this example, it would be dp/dn = 0. However, since I used did not use the PPE approach, I just imposed a condition of p=0 (or p=const since we are only interested in the gradient of p) at one of the boundaries (I believe the top-left point on the boundary).

Good day,

This is quite an interesting blog. I am studying for an M and right now I am about to complete the course work. I am now left with project work in order to get the M. So I have been thinking to do a Furnace Flow analysis, primarily looking at heat distribution. Now furnace consist of a domain B including the interior where the reactions take place and shell — being the furnace housing. I am tempted to thing think that part of the shell of course will form the domain boundary.

Now where you guys can be of help is the following, I have some experience in FEA and I comfortable with writing FEA codes (in Matlab), But the Question that remains is, Can I use FEA to solve the NS equations?

On the inside, I am worried about convection terms, on the walls the conduction lends itself some similarity with structural hookes law KU = F, so I suppose on the shell I am slightly covered.

So advice me on what to do.. I havent taken CFD courses however I have seen this problem being solved using CFD … I am hoping that the form of problem allows one to be able to tackle it with FEA.

Finite Element Analysis (FEA) is just one of the many numerical methods that can be used to solve complex problems. In the blog, I used finite different methods. However, the problem can indeed be solved using FEM. In the end, you will end up with a system of nonlinear algebraic equations that must be solved iteratively.

Can I get a copy of the code you show in this post? I’d have an immediate use for it, and it’s a bit long to re-type. (And it’s not possible to copy and paste from the website)

Thanks!

There is a link at the end of the blog that provides a link to downloading the CDF document of the blog. You can get the entire code from there.

Paritosh

I really enjoyed your post and now trying to extend it to annular flow case.I think it is better that you apply the zero velocity condition to bndryComp points rather than all boundary points and then overriding the top points horizontal velocity. You know, in the code you have defined bndryComp, but it has not been used.

Please explain the boundary condition setting part

I would try replacing the boundary conditions with something like the following:

inletVelocity[y_] := (*y*(1-y)*) Erf[10*y] + Erf[-10*(y - 1)] – 1;

eqnu[[top]] = uvar[[top]]; eqnu[[bot]] = uvar[[bot]];

eqnv[[top]] = vvar[[top]]; eqnv[[bot]] = vvar[[bot]];

eqnu[[right]] = dx[[right]].uvar; eqnv[[right]] = dx[[right]].vvar;

eqnu[[left]] = uvar[[left]] – inletVelocity[xygrid[[2]][[left]]];

eqnv[[left]] = vvar[[left]];

eqncont[[right]] = pvar[[right]];

the BCs above are: u(0,y)=inlet_profile. v(0,y)=0, du/dx(left,y) = dv/dx(left,y) = 0, u(top)=v(top)=0, u(bottom)=v(bottom)=0.

I would recommend replacing the BCs with one in the function DrivenCavitySolver and try for different Reynold’s number. Hope this helps!

say i have top and bottom boundary as no slip and inlet velocity = 0.001 ms-1 and outlet gauge pressure = 0

Sir

How to create a streamline and pathline for a rectangular closed cavity with one inlet gate and fill the rectangular cavity with viscous fluid.

Pls help me with some example to create pathline for filling the cavity in 10 sec and for every second, how the fluid flow inside cavity.

Regards

Rajashekaran K

For steady state flows, streamlines, streaklines, pathlines are all identical. If you want just streamlines, then I would recommend using the Mathematica functioin StreamPlot. For pathlines, if you are looking for animations, in the blog I have detailed how you can visualize the pathlines using NDSolve. I would recommend using that as a starting point for your analysis.

I was wondering what version of Mathematica you used. I am using Mathematica 9 and tried to follow this and got all the way down to where you use the FIndRoot function to solve the system. I received a message saying that the line search decreased the step size to within a tolerance specifed by AccuracyGoal and PrecisionGoal. My professor said it might have something to do with the version of Mathematica. If this is not the problem, do you have any suggestions to make this run better?

I might have used a different release of version 9 when I ran the code. FindRoot is a root-finding function and like all root-finding functions, it needs an initial condition. I would recommend playing around with the initial guess as a first step. For example, use an initial guess of -1 (instead of 1) and try running again.

sol = allVars /. FindRoot[govExpr, Thread[{allVars, -1}]] ;

I would like to point out that as the Reynold’s number increases, the initial guess can become quite sensitive. A proper initial guess in that sense is quite important.

The second thing to consider is grid convergence. Having a coarse grid can very well cause convergence issues. Also the difference order of the finite difference stencil. High difference orders could cause oscillations which prevent convergence. So, you can also play with these parameters and see the effect on solution and convergence.

I’ve run your code on MM v8, v9, and v10 at Re=1000, Difference Order = 2, and, since the solution was quicker, initial condition = -1. Version 8 and 9 took less than 2 minutes to calculate u, v, and p using FindRoot on a 100×100 grid. Version 10 took 75 s on a 20×20 grid, 262 s on 30×30, 3690 s on 40×40, and so many hours on 50×50 that I aborted it. Same problem for other Re and initial conditions and Difference Order. It seems to me that something must have been changed in FindRoot for v10. Can you confirm that and suggest how to deal with it? It’s a very nice algorithm but I can’t get it to complete for useful size grids on v10. I’ve kept v9 just for the purpose of running it. I’m running on a Mac, OS X 10.8.5, 2.66 GHz, 12 GB memory.

Dear Tom L

Re your post of Dec 22, 2014, I also find the exact same problem when moving from Mathematica Version 9 to Mathematic Version 10 . Under v9 grid sizes 100 by 100 took a few minutes, under v10 the calculation does not finish in hours. It’s a great program by Paritosh, I’d like to run it under v10.

Anthony Purcell

One approach that I can think of that will make FindRoot return a result in reasonable time is to generate initial conditions from the solution of a lower Reynolds number result. So, for example, to get the result at Reynold’s number 1000 with a 100×100 grid, one can do the following:

{u, v, p} = DrivenCavitySolver[100, {{0, 1}, {0, 1}}, {100, 100}, 4, "InitialGuess" -> 0];

xygrid = XandYGrid[{{0, 1}, {0, 1}}, {100, 100}];

grid = Flatten[Outer[List, Sequence @@ xygrid], 1];

ic = Join[Apply[u, grid, {1}], Apply[v, grid, {1}], Apply[p, grid, {1}]];

{u, v, p} = DrivenCavitySolver[1000, {{0, 1}, {0, 1}}, {100, 100}, 4, "InitialGuess" -> 1.5*ic]

If you want to get the result at 2000, then use the result from 1000, and so on. This process should allow the iterations to converge.

Hope this helps.

A reasonable approach but it didn’t solve the problem. v9 completed quickly as before but v10 took hours for a 50×50 so I aborted it. Something must have changed in FindRoot for v10 that severely degraded its performance. Perhaps the developer responsible for FindRoot could look into and resolve the problem for the next version. Thanks for your efforts.

Amazing example!

I am curious about a 3D case, say a cylinder. The air flows in through a small hole in the cylindrical surface and flows out through the two ends of the cylinder.

Jiri