Wolfram Computation Meets Knowledge

Using Mathematica to Simulate and Visualize Fluid Flow in a Box

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

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

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.

Comments

Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.

80 comments

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

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

      Reply
      • 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.

        Reply
      • Please post on 2d transient heat conduction

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

          Reply
          • While solving the transient heat conduction equation I get the following error:
            Warning: boundary and initial conditions are inconsistent.
            Please advice.

          • The message is not an error. It just means that the the initial condition values at the boundary don’t match the prescribed boundary condition values at time t = 0. As the simulation progresses, the result becomes consistent. NDSolve will still yield a valid result. I would recommend gonig through the NDSolve documentation to get a better idea of this behavior.

  2. 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!

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

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

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

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

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

    Reply
  8. 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

    Reply
  9. 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.

    Reply
    • 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.

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

        Reply
        • 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.

          Reply
  10. 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.

    Reply
    • 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.

      Reply
  11. Interesting post. What are the boundary conditions for P ?

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

      Reply
  12. 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.

    Reply
    • 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.

      Reply
  13. 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!

    Reply
  14. 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

    Reply
    • 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.

      Reply
  15. 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?

    Reply
    • 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.

      Reply
  16. 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.

    Reply
  17. 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

    Reply
    • 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.

      Reply
  18. 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.

    Reply
  19. 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

    Reply
  20. Good morning,

    Is there the possibility to extend the code for 3d flows in turbulent regime?
    If yes, how?

    Many thanks

    Reply
    • The construction of the matrices using finite-differences and the specification of the boundary conditions can be easily extended to 3D. I would recommend going through the advanced tutorial tutorial/NDSolveMethodOfLines in the documentation. Look at the section “Spatial derivative Approximation”. However, the formulation for solving the system will not hold. In the blog, we are solving a 2D steady-state flow. In 3D, turbulent regimes, we are looking at unsteady flows. Two very different problem. For the 3D turbulent case, the simplest method would be to make use of a fractional-step formulation. There is plenty of literature available online describing fractional-step methods.

      Reply
  21. Dear Tom L
    Thank you so much for all the information and explanation. I am a graduate student in civil engineering. I am still in my first step of coding and programming. However, I got valuable information and experience from what you wrote and explained. I tried to to compile your steps to run the code and visualize the flow in 2d but I got problems and errors. So. please if you can send me a copy from this code to find my error. Thanks again for you help.

    Reply
  22. I’m wondering if anyone of you had achieved to run this very nice program with Mathematica v10 in a reasonable amount of time, and, with which assumptions? Could it be Precision/AccuaryGoal/MaxIterations in the FindRoot?

    Best regards,
    Jürgen

    Reply
  23. Dear Paritosh, many thanks for your paper, very useful for my work on CFD (fans and combustion head flows). I would ask you: what about CFD simulations by finite differences and Reynolds number about 10^5 (such as air)? I have used your code or Picard iterations but for numerical stability and convergence the spatial steps for grid should be so small that unknowns became 10^5 or 10^6, a big problem for RAM and CPU. Do you have some ideas to suiggest? Thanks.

    Gianluca Argentini, RD& Dept., Riello Burners – Italy

    Reply
    • For such large Reynolds numbers, it is not surprising that you have 10^6 unknowns. Short of running such simulations on massive parallel clusters, one suggestion I can think of was to use domain decomposition methods to basically divide the domain up into small chunks and iterate/update between the various domains. There are several papers to that effect you can find in the literature. It will be slow but with such methods, at-least you will be able to control the RAM usage. Hope this helps.

      Reply
  24. Can you please let me know, how can I run this in Re = 20000 ??

    Reply
  25. Respected Paritosh Mokhasi,
    Hope you are doing well. I recently start working on this box flow. But i got some error while running the code in Mathematica 10 or 8. Can you suggest me where is the mistake. i got this error at

    eqnsU[[boundaries]] = varsU[[boundaries]]
    Part::partd: Part specification varsU[[{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,<>}]] is longer than depth of object. >>
    I really appreciate your effort and time. I can send you my code and my email address is ( maleeha.atlas29@gmail.com)

    Reply
  26. This simulation gives us several practical conditions. Here I inttended to highlighy only one. In the first square, we can imagina an open channel flow. We can observe the two inferior corners. We see some vortex there. This is the point where wh cut sharpened angles 90° an we put some pendant 45° instead. In this manner we improve the behaviour of the chammel and prlongue ots lifetime. Mathematics an numerical analysis, provide us a very nice comprehention of the phenomena and gives us the chance to invent devices to improve the behavioural of the structures an devices, mainly thinking in hydraulic facilities. Nice work. Congratulations! Transitories is one of the most difficult items on Hydraulics.

    Reply
  27. Sir
    Using the same method can we simulate space plasma behaviour

    Reply
  28. Thanks for this wonderful tutorial -it’s both extremely helpful and wonderfully informative. The biggest bottleneck for me was in deciphering the boundary conditions you applied. If I wanted to apply a Neumann/Robin boundary condition, what would be the best way to implement that?

    Reply
  29. I figured it out! I used compile to boost the calculation in the following manner :
    1) changed the line : {deqns, dvars} = {Join[eqnu, eqnv, eqncont], Join[uvar, vvar, pvar]};
    to
    dvars = Join[uvar, vvar, pvar];
    deqns = Compile[dvars, deqns,
    RuntimeOptions -> “EvaluateSymbolically” -> False];

    2) change Map[…] at the end by :
    Map[Compile[{x, y}, Interpolation@Join[grid, Transpose@List@#, 2] &,
    Partition[sol, Length[grid]],
    RuntimeOptions -> “EvaluateSymbolically” -> False]

    Now I solve a 100×100 in less than 6s against a time where I did not have the patience to wait for the result.

    for a 10×10 original 1.377304 s my version with compile 0.140221 s
    for a 20×20 original 110.167186 s my version with compile 0.308927s

    I hope that helps people who wants to use this code.

    BTW I am using Mathematica 10.3.0.0 on a MacBook Pro

    Reply
  30. interesting article. Want to try out soon.

    Reply
  31. sir can you post the code for generating pathlines? they are different than streamlines and vector field.

    Reply
  32. Respected Sir,
    Can I solve the following Non-linear partial differential equations directly.
    R = 0.5; M = 0.2;
    NDSolve[{D[f[x, y], x, x] + D[f[x, y], y, y] == -R * M* D[h[x, y], x],
    D[f[x, y], y]*D[h[x, y], x] – D[f[x, y], x]*D[h[x, y], y] == D[h[x, y], x, x] + D[h[x, y], y, y],
    f[0, y] == 0, h[0, y] == 0.5, f[1, y] == 0, h[1, y] == -0.5, f[x, 0] == 0, f[x,1]==0, h[x,0]==0,
    h[x,1]==0.

    Reply
    • These kind of elliptic nonlinear coupled PDEs cannot be solved directly within NDSolve currently. We hope to have this capability within NDSolve soon. In the meantime, you can use the steps provided in the blog to solve the system. Here is the code:

      pts = 100;
      nx = N@Range[0, 1, 1/(pts – 1)];
      grid = Flatten[Outer[List, nx, nx], 1];
      nL = Length[grid];
      psiv = Table[Unique[psi$], {i, nL}];
      thetav = Table[Unique[theta$], {i, nL}];

      {ddx, ddy, d2x, d2y} = NDSolve`FiniteDifferenceDerivative[#, {nx, nx}, “DifferenceOrder” -> 4][“DifferentiationMatrix”] & /@ {{1, 0}, {0, 1}, {2, 0}, {0, 2}};
      Lap = d2x + d2y;
      eye = IdentityMatrix[nL, SparseArray];

      {lbi, rbi, bbi, tbi} = Flatten[Position[grid, #]] & /@ {{0., y_}, {1., y_}, {x_, 0.}, {x_, 1.}};
      allbi = DeleteDuplicates[Join[lbi, rbi, bbi, tbi]];

      R = 0.5; M = 0.2;
      e1 = Lap.psiv + (R*M)*(ddx.thetav);
      e1[[allbi]] = eye[[allbi]].psiv;

      e2 = (ddy.psiv)*(ddx.thetav) – (ddx.psiv)*(ddy.thetav) – Lap.thetav;
      (*
      e2[[bbi]] = ddy[[bbi]].thetav;
      e2[[tbi]] = ddy[[tbi]].thetav;
      *)
      e2[[bbi]] = eye[[bbi]].thetav;
      e2[[tbi]] = eye[[tbi]].thetav;
      e2[[lbi]] = (eye[[lbi]].thetav) – 0.5;
      e2[[rbi]] = (eye[[rbi]].thetav) + 0.5;

      finaleqns = Join[e1, e2];
      finalvars = Join[psiv, thetav];

      sol = finalvars /. FindRoot[finaleqns, Thread[{finalvars, 0.5}]];

      Remove[“psi$*”, “theta$*”];
      {psiSol, thetaSol} = Map[Interpolation@Join[grid, Transpose@List@#, 2] &, Partition[sol, nL]];
      ContourPlot[#[x, y], {x, 0, 1}, {y, 0, 1}, Contours -> 20, ColorFunction -> “TemperatureMap”, ImageSize -> 400] & /@ {psiSol, thetaSol}

      If you want Neumann conditions on the top and bottom wall then uncomment the appropriate lines (and comment the two lines below it).

      Reply
  33. Respected Sir,
    Can I solve the following Non-linear partial differential equations directly.
    R = 0.5; M = 0.2;
    NDSolve [{D[f[x, y], x, x] + D[f[x, y], y, y] == -R * M* D[h[x, y], x],
    D[f[x, y], y]*D[h[x, y], x] – D[f[x, y], x]*D[h[x, y], y] == D[h[x, y], x, x] + D[h[x, y], y, y],
    f[0, y] == 0, h[0, y] == 0.5, f[1, y] == 0, h[1, y] == -0.5, f[x, 0] == 0, f[x,1]==0, h[x,0]==0,
    h[x,1]==0.
    The physical quantities of interest are the local Nusselt numbers [Nu]1 and [Nu]r , which are
    defined as
    [Nu]1=-kmnf/kf (∂θ/∂x)(x=0) , [Nu]r= -kmnf/kf (∂θ/∂x)(x=1).
    and the average Nusselt numbers [Nu]1 and [Nu]r, which are given by
    ([Nu]1 ) ̅=∫_0^1[Nu]1 dy , ([Nu]r ) ̅=∫_0^1[Nu]r dy.

    Reply
  34. Can you modify this code for a square lid-driven cavity with a circular hole at the centre?

    Reply
  35. This is a useful worked example.

    Do you know of any similar worked examples that include the time domain (I.e., t,x,y)? For example, the heat equation with additional non-linear terms, or spatially dependent diffusivities.

    I’m wondering about the feasibility of using the method of lines with FiniteDifferenceDerivative on a tensor grid. It would be nice to explore with a working example.

    Reply
    • Thanks for reaching out, Craig. There are several examples in the NDSolve/NDSolveValue documentation that deal with solving non-linear PDEs. Internally NDSolve does make use of method of lines with finite difference to integrate PDEs. Here is a example with nonlinear terms and spatial diffusive terms:

      sol = NDSolveValue[{D[u[x, t], t] + (1 + u[x, t])*D[u[x, t], x] == (Cos[x]^2) D[u[x, t], x, x],
      u[x, 0] == Sin[Pi x], u[0, t] == Sin[t], u[5, t] == 0}, u, {x, 0, 5}, {t, 0, 10}];
      Plot3D[sol[x, t], {x, 0, 5}, {t, 0, 10}, PlotRange -> All]

      I would highly recommend the visiting the documentation page tutorial/NDSolveOverview. There is an entire section describing the details of method of lines.

      Reply
  36. Hi there!
    My fourth year thesis project and I are working on modelling an irrotational vortex in Mathematica and was wondering if this model would be appropriate to work with. If you have any suggestions, please let me know!

    Reply
    • Hi Tess! Thanks for your question. Irrotational flow is in some sense simpler than solving the full-blown Navier-Stokes equations (which is what we are doing in this blog-post). Irrotational flow is characterized by flow-fields without vorticity, but the flow in a box result contains vorticity. You would need to look at potential flow theory for modelling irrotational vortex and perhaps look into vortex-panel methods if your project involves complex boundaries.

      Reply
  37. I’ve tried the above cdf in Mathematica v12, and it doesn’t seem to run in any reasonable time. It’s an interesting post. Would it be possible for you to modify/update the cdf so that it runs? I’m just starting to learn some hydrodynamics, so I’m afraid I’m unable to do so.

    I know you guys are busy, so no hard feelings if you can’t get to it. Thanks.

    Reply
  38. I also downloaded the post as cdf in Mathematica v11.3.0, and noticed that “In[58] ListAnimate[Table[Graphics[Thread…”, only seems to execute 1 step. On the next step one core goes to 100% and execution takes a long time. The only way to stop it is to hit the Pause button [||] in “Out[25]” animation control for points plot. Otherwise the statement never ends. if you terminate Mathematica task, it creates some kind of hang up or corruption that gets errors with “Evaluate Notebook”. Same thing happens at “Out[63]” for orange streamlines chart. What am I doing wrong?
    This does not seem to have anything to do with the link to “DrivenCavitySolver”.
    I will report later on my efforts to time those algorithms,

    Reply
  39. Can Mathematica handle the extension of this problem into three dimensions?

    Reply
    • Hi Virgilio,

      Mathematica can handle 3D case for flow in a box. The restriction will of-course be the amount of memory on your computer. Since the system of nonlinear equations becomes significantly more, the computation time needed to obtain the solution will also increase considerably.

      – Wolfram Blog Team

      Reply
  40. 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!

    Reply