Introduction

In addition to studying computational fluid dynamics (CFD), which is a study related to simulation of fluids such as air and water, I would like to summarize the knowledge necessary for constructing computational fluid dynamics (CFD) for water (in multiple articles).

I would like to write it so that even beginners can easily understand it. It seems that there are many mistakes, so please contact us if you find it. Also, I would appreciate it if you could comment on what is difficult to understand. We will update it from time to time.

Target audience

People who can use Python
Those who are interested in numerical calculation
Those who are interested in fluid mechanics
Those who understand basic university physics and mathematics (differential equations?)

series

Chapter 1: This article
Chapter 2: [Python] Fluid Simulation: Implementing Diffusion Equation
Supplementary position in Chapter 2: [Python] Article that enables high-speed sparse matrix calculation
Chapter 3: [Python] Fluid Simulation: From Linear to Nonlinear
Chapter 4: [Python] Fluid Simulation: Incompressible Navier-Stokes Equation

Rough content of this article

As a preliminary step to deal with the basic fluid equations required for water simulation, we will briefly summarize and implement the advection equations.

Change log

2020/02/19 There was an error in the text, so I corrected it. Thank you to @picric_acid for commenting!
2020/02/21 There was a bug in the code, so it has been fixed. Thank you to @fluxCapacitor for commenting!

Advection equation (equation that expresses the movement of matter)

We will solve the one-dimensional advection equation by the finite difference method.

Advection equation $ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $

This equation means that a function u moves at a velocity c. In short, it is an equation that expresses that a physical quantity moves along a flow called velocity c.

Discrete this and solve the equation. This time, we will solve it by a method called Positive Solution. The explicit method is simply a method of finding future values using only the current time value. In addition to the explicit method, there is also a method called the implicit method, which is a method of calculating based on the current value and the future value.

The explicit method to be implemented this time is as follows. I chose the following four at will.

FTCS(Forward in Time and Central Difference in Space)
Upwind differencing
Lax-Wendroff scheme
CIP(Constrained Interpolation Profile scheme) method

For each, we will write down the discretized expressions. In the expression, the subscript j means place and n means time. Therefore, $ f_j ^ n $ means the value of the function f located at the time $ n $ at the location $ x_j $.

FTCS(Forward in Time and Central Difference in Space) 2. $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) $
Upwind differencing 3. $ u_j^{n+1} = u_j^{n} - c \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j}^n - u_{j-1}^n \right) $
Lax-Wendroff scheme 4. $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) + \frac{c^2}{2} \left( \frac{\Delta t}{\Delta x} \right)^2 \left(u_{j+1}^n - 2 u_j^n+ u_{j-1}^n \right) $
CIP(Constrained Interpolation Profile scheme) method

Details will be described later. Alternatively, please refer to Introduction to CIP Law.

Stable conditions of explicit method

I would like to briefly mention the important stability conditions when using the explicit method. Numerically, the information propagation velocity $ \ frac {\ Delta t} {\ Delta x} $ (the velocity obtained by dividing the grid width by time) must be greater than or equal to the physical signal velocity $ c . Otherwise, the physical phenomenon will move faster than the numerical calculation, and the physical phenomenon cannot be expressed numerically. Based on that idea $ c \leq \frac{\Delta x}{\Delta t} $ The condition is derived. Transformed this $ \nu = c \frac{\Delta t}{\Delta x} \leq 1 $$ Is called the CFL condition and is used as a condition that is numerically stable. The left side is also called the Courant number. In short, the time step width $ \ Delta t $ should be less than or equal to $ \ frac {\ Delta x} {c} $. Note that these CFL conditions have no relation when using the implicit method, so the implicit method has the advantage that it can be made larger in time increments than the explicit method. However, the calculation becomes complicated accordingly.

Implementation

The language used is Python. I will write it without using many Numpy functions so that the calculation formula is easy to understand.

Calculation target

Calculate the advection of a square wave. The grid width $ \ Delta x $, the physical signal velocity $ c $, and the time step width $ \ Delta t $ are all calculated as constants. This time, it is calculated as $ \ Delta x = 1 $, $ c = 1 $, and $ \ Delta t = 0.2 $ according to the references. Therefore, it is fixed at $ CFL = 0.2 $ and the CFL condition is satisfied.

# Creat square wave
Num_stencil_x = 101
x_array = np.arange(Num_stencil_x)
u_array = np.where(((x_array >= 30) |(x_array < 10)), 0.0, 1.0)
u_lower_boundary = 0.0
u_upper_boundary = 0.0
Time_step = 200
Delta_x = max(x_array) / (Num_stencil_x-1)
C = 1
Delta_t = 0.2
CFL = C * Delta_t / Delta_x
total_movement = C * Delta_t * (Time_step+1)
exact_u_array = np.where(((x_array >= 30 + total_movement) |(x_array < 10 + total_movement)), 0.0, 1.0)
plt.plot(x_array, u_array, label="Initial condition")
plt.plot(x_array, exact_u_array, label="Answer")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

FTCS(Forward in Time and Central Difference in Space) $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) $

The solution vibrates. I don't get the answer at all

u_ftcs = u_array.copy()
#Set time step
for n in range(Time_step):
    u_old = u_ftcs.copy()
    u_ftcs[0] = u_old[0] - CFL / 2 * (u_old[1] - u_lower_boundary)
    u_ftcs[-1] = u_old[-1] - CFL / 2 * (u_upper_boundary - u_old[-1])
    for j in range(1,Num_stencil_x-1, 1):
        u_ftcs[j] = u_old[j] - CFL / 2 * (u_old[j+1] - u_old[j-1])
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_ftcs, label="FTCS")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(FTCS)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

2. Upwind differencing

u_j^{n+1} = u_j^{n} - c \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j}^n - u_{j-1}^n \right)

The numerical diffusion is large and the solution becomes smooth.

u_upwind = u_array.copy()
#Set time step
for n in range(Time_step):
    u_old = u_upwind.copy()
    u_upwind[0:1] = u_old[0] - CFL * (u_old[1] - u_lower_boundary)
    for j in range(1,Num_stencil_x):
        u_upwind[j:j+1] = u_old[j] - CFL * (u_old[j] - u_old[j-1])
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_upwind, label="Upwind differencing")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(Upwind differencing)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

Lax-Wendroff scheme $ u_j^{n+1} = u_j^{n} - \frac{c}{2} \left( \frac{\Delta t}{\Delta x} \right) \left(u_{j+1}^n - u_{j-1}^n \right) + \frac{c^2}{2} \left( \frac{\Delta t}{\Delta x} \right)^2 \left(u_{j+1}^n - 2 u_j^n+ u_{j-1}^n \right) $

There is numerical vibration and blunt solution.

u_lw= u_array.copy()
#Set time step
for n in range(Time_step):
    u_old = u_lw.copy()
    u_lw[0:1] = u_old[0] - CFL / 2 * (u_old[1] - u_lower_boundary) \
                    + CFL**2 / 2 * (u_old[1] - 2 * u_old[0] + u_lower_boundary)
    u_lw[-1] = u_old[-1] - CFL / 2 * (u_upper_boundary - u_old[-1]) \
                    + CFL**2 / 2 * (u_upper_boundary - 2 * u_old[-1] + u_old[-2])
    for j in range(1,Num_stencil_x-1, 1):
        u_lw[j] = u_old[j] - CFL / 2 * (u_old[j+1] - u_old[j-1]) \
                    + CFL**2 / 2 * (u_old[j+1] - 2 * u_old[j] + u_old[j-1])
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_lw, label="Lax-Wendroff scheme")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(Lax-Wendroff scheme)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

CIP(Constrained Interpolation Profile scheme) method

The basic idea of the CIP method is that when a function f is discretized, the physical quantities of a point $ x_j $ and a point $ x_ {j-1} $ next to it are smoothly connected. By doing this, the differential values of the respective physical quantities $ f_j $ and $ f_ {j-1} $ can be used for time evolution without diverging the slope between the lattices.

In other words, under the condition that the ** function value and the differential value are continuous on the grid points **, the relationship between the two grid points can be expressed by the cubic complementary function $ F (x) $.

In the area of $ x_ {j-1} <= x <x_j $ $ F(x) = a X^3 + b X^2 +c X + d \quad ,where \quad X = x - x_j $ If $ F (x) $ is defined by cubic completion as in, ** from the condition that the function value and the differential value are continuous on the grid point **

F(0) = f_j^n, \quad F(\Delta x) = f_{j+1}^n, \quad \frac{\partial F(0)}{\partial x} = \frac{\partial f_j^n}{\partial x}, \quad \frac{\partial F(\Delta x)}{\partial x} = \frac{\partial f_{j+1}^n}{\partial x} \\ , where \quad \Delta x = x_{j+1} - x_j

Is satisfied, so if you substitute this into the formula of $ F (x) $ and find $ a, b, c, d $,

a = \frac{\partial f_j^n / \partial x + \partial f_{j+1}^n / \partial x }{\Delta x^2} + \frac{2 (f_j^n - f_{j+1}^n)}{\Delta x^3} \\
b =  \frac{3 (f_{j+1}^n - f_j^n)}{\Delta x^2} - \frac{2( \partial f_j^n / \partial x + \partial f_{j+1}^n / \partial x )}{\Delta x} \\
c = \frac{\partial f_j^n}{\partial x}\\
d = f_j^n

From the above, when the velocity c is constant, $ \ frac {\ partial u} {\ partial t} + c \ frac {\ partial u} {\ partial x} = 0 $ also satisfies the differential value, so the time step width $ \ Delta t (n + 1 step) The value and derivative value after $ seconds are

u_j^{n+1} = F(x_j - c \Delta t)= a \xi^3 + b \xi^2 + c\xi + d \\
\frac{\partial u_j^{n+1}}{\partial x} = \frac{d F(x_j - c \Delta t)}{d x} = 3 a \xi^2 + 2 b \xi + c \\
, where \quad \xi = - c \Delta t

Numerical diffusion and vibration are overwhelmingly smaller than the other three methods. ~~ However, the advection was not as good as References. ~~ (2020/02/20 After making the corrections I received in the comments, it has improved considerably.)

u_cip= u_array.copy()
partial_u_cip = ((np.append(u_cip[1:], u_upper_boundary) + u_cip)/2 - (np.append(u_lower_boundary, u_cip[:-1]) + u_cip)/2)/ Delta_x
#Evolve time
for n in range(Time_step):
    u_old = u_cip.copy()
    partial_u_old = partial_u_cip.copy()
    u_cip[0] = 0
    partial_u_cip[0] = 0
    for j in range(1, Num_stencil_x):
        a = (partial_u_old[j] + partial_u_old[j-1]) / Delta_x**2 - 2.0 * (u_old[j] - u_old[j-1]) / Delta_x**3
        b = 3 * (u_old[j-1] - u_cip[j]) / Delta_x**2 + (2.0*partial_u_old[j] + partial_u_old[j-1]) / Delta_x
        c = partial_u_old[j]
        d = u_old[j]
        xi = - C * Delta_t  # C > 0
        u_cip[j:j+1] = a * xi**3 + b * xi**2 + c * xi + d
        partial_u_cip[j:j+1] = 3 * a * xi**2 + 2 * b * xi + c
plt.plot(x_array, exact_u_array, label="Answer")
plt.plot(x_array, u_cip, label="CIP method")
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("u(CIP)")
plt.xlim(0, max(x_array))
plt.ylim(-0.5,1.5)

Summary

The basics of the advection equations that make up the fluid equation are summarized with a brief implementation. The following four have been implemented.

FTCS(Forward in Time and Central Difference in Space)
Upwind differencing
Lax-Wendroff scheme
CIP(Constrained Interpolation Profile scheme) method

Next time, I will talk about Diffusion equation.

References

http://zellij.hatenablog.com/entry/20140805/p1
https://www.cradle.co.jp/media/column/a233
http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%B0%DC%CE%AE%CB%A1
http://www.astro.phys.s.chiba-u.ac.jp/netlab/mhd_summer2001/cip-intro.pdf
http://www.astro.phys.s.chiba-u.ac.jp/netlab/summer-school_2004/TEXTBOOK/text3-1.pdf

[Python] Fluid simulation: Implement the advection equation