Jorge Cardona's Blog Mathematics

[DRAFT] Computational Fluid Dynamics

Updated: December 22, 2024

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 u of a PDE, is usually a function defined on a space-time domain [0,T]×Ω. Boundary conditions are conditions that must be satisfied by u on the boundary of Ω. The simplest situation is to force the function to vanish, i.e. u(t,x)=0 for every xΩ and every t[0,T].

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 ut=2ux2 models the transfer of heat on a given domain. We will restrict ourselves to a one-dimensional rod represented by the domain [0,1].

To construct our simulation we will consider the initial condition given by u(0,x)=12πσ2ex2/2σ2.

The periodic boundary conditions simplify the problem of approximating the second spatial derivative.

A solution u:[0,T]×[0,1]R 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 Δt in time and Δx in space. Then, using the numpy notation to access the elements of a tensor the approximation is u[i,j]=u(iΔt,jΔx).

In sake of simplicity we will consider periodic boundary conditions. In this case the second derivative is defined via a central differences as uxx(iΔt,jΔx)=u[i,j1]+u[i,j+1]2u[i,j]Δx2 we can implement this in python as follows, the simplest simulation is just an approximation using Euler method as follows:

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
import numpy as np

def solve_heat(U_0, dT, dX, steps=10, skip=0):

# If this coefficient is small the simulation is stable.
r = dT / dX ** 2

# Initialize first the variables.
U = np.zeros((steps, X.shape[0]))
T = np.zeros(steps)

# Initial conditions.
U[0] = U_0

# Start the iterations.
for i in range(1, steps):

U_ = U[i - 1]

for j in range(skip + 1):
U_ = (1 - 2 * r) * U_ + r * (np.roll(U_, 1) + np.roll(U_, -1))

U[i] = U_
T[i] = T[i -1] + dT * (skip + 1)

return T, U

def initial_condition(x, sigma=0.2):
return np.exp(-x**2 / (2 * sigma**2)) / np.sqrt(2 * np.pi * sigma ** 2)

# Spatial domain.
X = np.linspace(-1, 1, 123)

# Initial condition.
U_0 = initial_condition(X)

# Solve the equation.
T, U = solve_heat(U_0, dT=0.0001, dX=X[1] - X[0], steps=40, skip=100)

# Generate a GIF.
generate_gif(X, U, T, 'evolution_heat_1d_2.gif')

gif

Burgers equation in one dimension.

The Burgers equation is the following equation ut+uux=0

If we try to simulate this in the naive way

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
def solve_burgers(U_0, dT, dX, steps=10, skip=0):

# If this coefficient is small the simulation is stable.
r = dT / dX

# Initialize first the variables.
U = np.zeros((steps, X.shape[0]))
T = np.zeros(steps)

# Initial conditions.
U[0] = U_0

# Start the iterations.
for i in range(1, steps):

U_ = U[i - 1]

for j in range(skip + 1):
U_ = U_ * (1 - r * (np.roll(U_, -1) - U_))

U[i] = U_
T[i] = T[i -1] + dT * (skip + 1)

return T, U

def initial_condition(x, sigma=0.2):
return np.exp(-x**2 / (2 * sigma**2)) / np.sqrt(2 * np.pi * sigma ** 2)

# Spatial domain.
X = np.linspace(-1, 1, 123)

# Initial condition.
U_0 = initial_condition(X)

#
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')

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

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.

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
def solve_burgers(U_0, dT, dX, mu=0.0, steps=10, skip=0):

# If these coefficient are too large the simulation is not stable.
r = dT / dX
r2 = dT / dX ** 2

# Initialize first the variables.
U = np.zeros((steps, X.shape[0]))
T = np.zeros(steps)

# Initial conditions.
U[0] = U_0

# Start the iterations.
for i in range(1, steps):

U_ = U[i - 1]

for j in range(skip + 1):
U_ = U_ * (1 - r * (np.roll(U_, -1) - U_) - 2 * r2 * mu) \
+ mu * r2 * (np.roll(U_, 1) + np.roll(U_, -1))

U[i] = U_
T[i] = T[i -1] + dT * (skip + 1)

return T, U

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

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.

Two dimensional problems

So, we have seen that simpl e