Wolfram Blog
Paritosh Mokhasi

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.

da Vinci sketch

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:

Motion of a fluid

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/m2) represents the pressure; ν (m2/s) is the viscosity of the fluid; and ρ (kg/m3) 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.

Flow in a box sketch

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:

Navier-Stokes equations simplified

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

Three unknowns and three equations

where a characteristic velocity U0(m/s) and a characteristic length L (m) are used, the equations become

Equations with velocity and length

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 xi, and the solution is to be found at u(xi).

2) The derivatives are expressed as some linear combination of u(xi).

3) The final system is changed from a system of differential equations to a system of algebraic equations for u(xi). 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:

Define a set of grid points

ListPlot of the PlotRange

At each of these grid points, the velocities u(xi, yi) and v(xi, yi) and pressure p(xi, yi) have to be computed. This means that each point contains three unknowns. We can give each unknown a name.

nL = Length[grid]; varsU = Table[Unique[u$], {nL}]; varsV = Table[Unique[v$], {nL}]; varsP = Table[Unique[p$], {nL}];

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

du/dx(Subscript[x, i],Subscript[y, i]) \[TildeTilde] \[Sum]Subscript[a, ij]u(Subscript[x, j],Subscript[y, j]).

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:

ugrid = h Range[0, 8]; Thread[{Thread[     Map[f'[#] &, ugrid] ==       NDSolve`FiniteDifferenceDerivative[Derivative[1], ugrid,        Map[f, ugrid]]]}] // Grid

ugrid

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).

first order derivative for Sin(x) where x \[Element](0,1)

{-0.0000196916,   4.90987*10^-6, -3.263*10^-6, -3.18067*10^-6, -3.06655*10^-6, \ -2.92179*10^-6, -2.74785*10^-6, -2.54644*10^-6, -2.31959*10^-6,   3.41884*10^-6, -0.0000134268}

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.

DifferenceOrder

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

Nonlinear algebraic equations

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:

Find the grid points that lie on the boundary AND the equation numbers associated with them

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

u(0,y)=u(1,y)=u(x,0)=0,  u(x,1)-1=0, v(0,y)=v(1,y)=v(x,0)=v(x,1)=0

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.

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 3n 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.

govExpr = Join[eqnsU, eqnsV, eqnsCont] ;  allVars = Join[varsU, varsV, varsP];  Length /@ {govExpr, allVars}

{30000, 30000}

The function FindRoot can be used to solve this system.

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

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.

A black-box continuous function

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

>LineIntegralConvolutionPlot

Very attractive flow visualization

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

A number of helper/utility functions in order to put the code in a resuable form

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

The main solver function

The main solver function cont.

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.

What happens to the flow when the Reynolds number is increased?

GraphicsGrid

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.

{u, v, p} = DrivenCavitySolver[2000, {{0, 1}, {0, 1}}, {100, 100}, 4];

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.

ipt = Outer[List, Range[0.05, 0.95, 0.9/19],     Range[0.05, 0.95, 0.9/19]]; ipt = Transpose[Flatten[ipt, 1]];

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.

A time window of 10 time-units

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

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.

Equation to release a continuous stream of dye at one location in space for a short period of time

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.

Change the aspect ratio

The vortex at the bottom-left of the box

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

The box is now made wide instead of narrow

Wide visualization

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.

Leave a Comment

26 Comments


Math Nerd

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.

Posted by Math Nerd    July 9, 2013 at 2:56 pm
    Paritosh Mokhasi

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

    Posted by Paritosh Mokhasi    July 10, 2013 at 9:13 am
      Asim

      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.

      Posted by Asim    July 10, 2013 at 2:48 pm
      Debashish Deka

      Please post on 2d transient heat conduction

      Posted by Debashish Deka    March 10, 2014 at 10:14 pm
        Paritosh Mokhasi

        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}]

        Posted by Paritosh Mokhasi    April 16, 2014 at 10:02 am
Nicholas Lucas

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!

Posted by Nicholas Lucas    July 9, 2013 at 3:30 pm
    Paritosh Mokhasi

    Thanks for catching the typo. I corrected the equation.

    Posted by Paritosh Mokhasi    July 10, 2013 at 8:55 am
Paco

Amazing!

Posted by Paco    July 9, 2013 at 5:14 pm
P. Fonseca

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.

Posted by P. Fonseca    July 9, 2013 at 11:02 pm
Maximilien Levesque

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

Posted by Maximilien Levesque    July 10, 2013 at 10:44 am
John Snyder

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

Posted by John Snyder    July 10, 2013 at 1:48 pm
Paritosh Mokhasi

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

Posted by Paritosh Mokhasi    July 12, 2013 at 10:28 am
Kevin

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.

Posted by Kevin    August 17, 2013 at 12:47 pm
syed osama ali

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

Posted by syed osama ali    March 8, 2014 at 11:44 am
Kirolosse GIRGIS

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.

Posted by Kirolosse GIRGIS    April 11, 2014 at 2:12 pm
    Paritosh Mokhasi

    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.

    Posted by Paritosh Mokhasi    April 16, 2014 at 10:09 am
      Kirolosse GIRGIS

      Isn’t the difficulty only how to impose the boundary conditions .. ?
      and how can I solve the problem of the irregular geometry..

      Posted by Kirolosse GIRGIS    April 23, 2014 at 1:50 am
        Paritosh Mokhasi

        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.

        Posted by Paritosh Mokhasi    April 30, 2014 at 3:32 pm
JD Metzger

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.

Posted by JD Metzger    August 11, 2014 at 9:41 am
    Paritosh Mokhasi

    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.

    Posted by Paritosh Mokhasi    August 13, 2014 at 10:28 am
Ivan

Interesting post. What are the boundary conditions for P ?

Posted by Ivan    September 10, 2014 at 8:21 pm
    Paritosh Mokhasi

    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).

    Posted by Paritosh Mokhasi    September 11, 2014 at 12:10 pm
ernest

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.

Posted by ernest    October 3, 2014 at 5:33 am
    Paritosh Mokhasi

    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.

    Posted by Paritosh Mokhasi    October 8, 2014 at 10:14 am


Leave a comment

Loading...

Or continue as a guest (your comment will be held for moderation):