# 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/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.

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 U0(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 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:

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.

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 NDSolveFiniteDifferenceDerivative that automatically computes these coefficients. Here is an example:

The function NDSolveFiniteDifferenceDerivative 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 NDSolveFiniteDifferenceDerivative, 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 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.

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.

 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
 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
 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
 Please post on 2d transient heat conduction Posted by Debashish Deka    March 10, 2014 at 10:14 pm
 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
 While solving the transient heat conduction equation I get the following error: Warning: boundary and initial conditions are inconsistent. Please advice. Posted by Vikalp Luthra    August 25, 2016 at 10:49 am
 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. Posted by Paritosh Mokhasi    September 6, 2016 at 10:12 am
 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
 Thanks for catching the typo. I corrected the equation. Posted by Paritosh Mokhasi    July 10, 2013 at 8:55 am
 Amazing! Posted by Paco    July 9, 2013 at 5:14 pm
 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
 I also never post, but this one is amazing. Thanks!!!! Posted by Maximilien Levesque    July 10, 2013 at 10:44 am
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 Interesting post. What are the boundary conditions for P ? Posted by Ivan    September 10, 2014 at 8:21 pm
 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
 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
 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
 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! Posted by Chris Reed    October 14, 2014 at 8:38 pm
 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. Posted by Paritosh Mokhasi    October 20, 2014 at 3:43 pm
 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! Posted by Paritosh Mokhasi    November 19, 2014 at 12:33 pm
 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 Posted by Raj    November 19, 2014 at 12:00 am
 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. Posted by Paritosh Mokhasi    November 19, 2014 at 3:01 pm
 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? Posted by Michael    November 27, 2014 at 1:38 pm
 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. Posted by Paritosh Mokhasi    December 1, 2014 at 10:37 am
 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. Posted by Tom L    December 22, 2014 at 12:01 pm
 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 Posted by Anthony Purcell    January 13, 2015 at 3:00 pm
 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. Posted by Paritosh Mokhasi    January 14, 2015 at 1:10 pm
 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. Posted by Tom L    January 18, 2015 at 11:24 pm
 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 Posted by Jiri R.    February 25, 2015 at 7:52 am
 Good morning, Is there the possibility to extend the code for 3d flows in turbulent regime? If yes, how? Many thanks Posted by Matteo    May 5, 2015 at 4:03 am
 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. Posted by Paritosh Mokhasi    May 6, 2015 at 10:07 am
 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. Posted by Zeyad Sulaiman    May 17, 2015 at 10:13 pm
 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 Posted by Jürgen    July 3, 2015 at 8:14 am
 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 Posted by Gianluca Argentini    September 17, 2015 at 6:39 am
 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. Posted by Paritosh Mokhasi    September 24, 2015 at 9:41 am
 Can you please let me know, how can I run this in Re = 20000 ?? Posted by Selva    October 8, 2015 at 5:48 pm
 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) Posted by Maleeha Atlas    October 13, 2015 at 4:49 am
 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. Posted by Dr. Eng. Carlos Rivera    October 21, 2015 at 11:54 pm
 Sir Using the same method can we simulate space plasma behaviour Posted by sijo    December 31, 2015 at 4:54 am
 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? Posted by Joe    January 30, 2016 at 2:37 pm
 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 Posted by Remington    February 29, 2016 at 11:57 am
 I am confused about your first use of Compile for deqns. When I make the changes you recommend the code does not evaluate and throws errors during the findroot stage. Can you please elaborate a little (or check if there are any typo’s in your original comment) thanks! Posted by Sander    May 19, 2016 at 11:09 pm
 interesting article. Want to try out soon. Posted by nirmal    May 3, 2016 at 6:22 am
 sir can you post the code for generating pathlines? they are different than streamlines and vector field. Posted by sadaf    November 17, 2016 at 11:40 pm
 For steady flows, streamlines, streaklines and pathlines are identical. So, the streamline plots would suffice. Vector field can be computed using the function VectorPlot. Posted by Paritosh Mokhasi    November 18, 2016 at 4:00 pm
 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. Posted by P. Sudarsana Reddy    October 28, 2017 at 8:15 am
 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} = NDSolveFiniteDifferenceDerivative[#, {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). Posted by Paritosh Mokhasi    October 31, 2017 at 9:35 am
 very intresting Posted by hanen Oueslati    November 28, 2017 at 3:00 am
 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. Posted by P. Sudarsana Reddy    August 23, 2018 at 3:55 am
 Can you modify this code for a square lid-driven cavity with a circular hole at the centre? Posted by Mr. Saha    September 17, 2018 at 12:00 am
 Good question! The code can be modified to have an arbitrary body inside the box by using the “Immersed Boundary Method”. However, the modifications will not be trivial (or at-least not one-liners). This will take a little bit of time to make the modifications. Posted by Wolfram Blog    September 18, 2018 at 9:34 am
 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. Posted by W Craig Carter    October 6, 2018 at 9:35 am
 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. Posted by Wolfram Blog    October 9, 2018 at 9:49 am