In this post I will present how to use python as a tool to simulate fluid mechanics. It will be very basic at first and I will try to increase the complexity either by continue in other posts or simply by expanding the current one.
The code in this post is automatically collected at generation time and run, to ensure that the code actually run and generate the correct plot. You can find the collected code here.
A word on visualization
It is not simple to visualize the solutions of a PDE, specially in three-dimensional problems. For one-dimensional problems a simple video of a the plot is enough, a video of a two-dimensional level-based coloring is good as well, for three-dimensional problems, we could just consider sections of the domain.
Since we are interested here in learning how to construct simulations of fluids we need to understand what it means to simulate a solution of a PDE. Without goint into much details, we are constructing approximations to the original problem, in particular we define a new equation in differences that approximate, in some sense, the original PDE.
We need to be careful with the approximation of the different partial derivatives appearing in the original equation and of the boundary conditions.
Boundary conditions
A solution of a PDE, is usually a function defined on a space-time domain . Boundary conditions are conditions that must be satisfied by on the boundary of . The simplest situation is to force the function to vanish, i.e. for every and every .
Boundary conditions capture different physical aspects of a problem in consideration, as well as changing the analysis required to find interesting results associated with a particular problem. In simulating a solution of a PDE there are also some practical considerations that change according to the boundary conditions.
Periodo
One-dimensional examples
The heat equation with periodic boundary conditions
The heat equation models the transfer of heat on a given domain. We will restrict ourselves to a one-dimensional rod represented by the domain .
To construct our simulation we will consider the initial condition given by
The periodic boundary conditions simplify the problem of approximating the second spatial derivative.
A solution of the heat equation is approximated by a two dimensional numpy tensor in python. We use a simple rectangular grid defined by the step sizes in time and in space. Then, using the numpy notation to access the elements of a tensor the approximation is .
In sake of simplicity we will consider periodic boundary conditions. In this case the second derivative is defined via a central differences as we can implement this in python as follows, the simplest simulation is just an approximation using Euler method as follows:
# T, U = solve_burgers(U_0, dT=0.0001, dX=X[1] - X[0], steps=20, skip=100) generate_gif(X, U, T, 'evolution_burgers_1d.gif')
This method is unstable for the Burgers equation, even if we decrease the time and space steps:
1 2 3 4
X = np.linspace(-1, 1, 1023) u_0 = initial_condition(X) T, U = solve_burgers(u_0, dT=0.0000001, dX=X[1] - X[0], steps=20, skip=10000) generate_gif(X, U, T, 'evolution_burgers_1d_2.gif')
Solutions to the burgers equation loose they intial regularity, and any discontinuity appearing in the simulation would be amplified rapidly. In comparison, adding a viscosity term the system comes back to be stable.
X = np.linspace(-1, 1, 123) u_0 = initial_condition(X) T, U = solve_burgers(u_0, dT=0.00001, mu=0.025, dX=X[1] - X[0], steps=80, skip=4000) generate_gif(X, U, T, 'evolution_burgers_vis_1d.gif')
The theory behind the Burgers’ equation is very interesting, but I will not go into it now, there are weak solutions which capture in a better way the problems and allow for a better understanding.