[Python] Fluid simulation: Incompressible Navier-Stokes equation

Introduction

In addition to studying computational fluid dynamics (CFD), I would like to summarize the knowledge necessary for constructing computational fluid dynamics (CFD) for liquids. It seems that there are many mistakes, so please contact us if you find it.

Target audience

series

Rough content of this article

This article deals with single-phase numerical simulations such as gas only or liquid only. After briefly explaining the governing equation of the fluid, it is configured to implement the numerical calculation algorithm of the incompressible fluid required for the simulation of the liquid. Finally, we will simulate the cavity flow (feeling of spinning the fluid in the area / figure below). Regarding the cavity flow, @ eigs's Numerical solution of incompressible Navier-Stokes equation 1: Introduction is easy to understand and the graph is also beautiful.

2d-cavity.gif

table of contents

chapter title Remarks
1. Fluid solution
1.1. Fluid governing equation
1.2. Compressible and incompressible
2. Incompressible fluid calculation algorithm
2.1. MAC method Marker And Cell method
2.2. Fractional Step method Also called partial step method
2.3. SIMPLE method Semi-Implicit Method for Pressure-Linked Equation method
2.4. CCUP method CIP-Combined Unified Procedure method
2.5. Spatial difference Regarding variable placement
3. Implementation CCUP method

1. Fluid solution

1-1. Fluid governing equation

A fluid consists of three equations: the continuity equation, the Navier Stokes equation, and the energy equation. Each equation is derived from mass conservation, momentum conservation, and energy conservation. It may be more confusing, but I'll also add an image of each formula.

a. Continuity equation

\frac{\partial \rho}{\partial t} +  (\vec{u} \cdot \nabla) \rho  =   - \rho \nabla \cdot \vec{u}  \\
\quad or \\
\frac{\partial \rho}{\partial t} + \nabla \rho \vec{u} = 0

It can be derived from the fact that the mass change within a certain control volume is equal to the surface integral of the mass flow rate $ \ rho \ vec {u} $ flowing in and out of the surface. To explain with an easy-to-understand example, the above formula shows that the amount of change in money in the wallet is equal to the amount obtained by subtracting the amount of money that has passed through the doorway of the wallet.

money.png

b. Navier Stokes equation (equation of motion of fluid)

\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + (\mu + \lambda) \nabla (\nabla \cdot \vec{u}) + \vec{F_B} \\
 \quad or \\
\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \nabla \cdot {\bf T}_{\nu} + \vec{F_B} \\
, where \quad {\bf T}_{\nu} = \lambda (\nabla \cdot \vec{u}) {\bf I} + 2 \mu {\bf D} \quad :Viscous stress tensor\\
{\bf D}  = \frac{1}{2} \left\{(\nabla \vec{u}) + (\nabla \vec{u})^T \right\} :Strain rate tensor\\
\quad \vec{F_B} \quad :Body force(Mainly gravity) \\
\quad \mu :Viscosity coefficient

The momentum of the fluid is $ pressure gradient: \ nabla p $, $ diffusion due to viscosity: \ nabla \ cdot {\ bf T} _ {\ nu} $, $ body force (gravity etc.): \ vec {F_B} $ Indicates that you will be affected. To give an image example, it may be helpful to imagine that when a boy kicks a ball, the speed of the ball changes due to the effects of pressure difference, viscosity, and gravity.

navier_stokes.png

c. Energy equation

 \rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = -p \nabla \cdot \vec{u} - \nabla \cdot (\kappa \nabla T) + {\bf T}_{\nu} : {\bf D} + \dot{L} \\
, where \quad  \dot{L} :Amount of heat generated by chemical reaction\\
\quad \kappa :Thermal conductivity

However,

{\bf T}_{\nu} : {\bf D} = tr({\bf T}_{\nu} \cdot {\bf D})

The amount of change in the internal energy of the fluid is the work $ \ nabla \ cdot \ vec {u} $ due to pressure, the energy balance due to heat conduction $ \ nabla \ cdot (\ kappa \ nabla T) $, and the dissipative energy due to viscosity $ {\ bf. It means that it is affected by $ \ dot {L} $ such as T} _ {\ nu} \ cdot \ nabla \ vec {u} $ and chemical reactions. This formula allows the temperature change of the fluid to be calculated. In short, the temperature change of a fluid is equal to the amount of heat applied.

energy_equation.png

Summary

The figure shown above can be summarized as follows.

flow_equations.png

The above is the governing equation of fluid in a non-conservative system, and if the right side is set to 0, it means the convection equation of density, velocity, and internal energy (an equation that means mass transfer). A non-conservative system is a form in which variables exist outside the derivative, and there is no guarantee that the calculation results satisfy the conservation law, so it is necessary to confirm the conservativeness when performing numerical calculations. The opposite of the non-conservative system is called the conservative system, which is an expression in which all variables are included in the derivative.

1-2. Compressible and incompressible

The equation shown in Section 1-1 is the governing equation for compressible fluids. A compressible fluid is a fluid whose flow velocity is about 0.3 times the speed of sound (Mach number is 0.3), and a fluid with a slower flow is generally recognized as an incompressible fluid. In short, if the fluid moves too fast, the fluid itself expands and contracts (changes in density), making it compressible, and if it is slow or dense, the fluid does not expand or contract and is treated as incompressible. Therefore, the liquid is basically treated as uncompressed.

flow.png

The formula shown in Section 1-1 is transformed into an uncompressed formula. All we have to do is treat the density as a constant rather than a variable to simplify the equation. The Navier Stokes equation and the energy equation are simplified by using the continuity equation. In addition, we will ignore the effects of chemical reactions here.

a. Continuity equation

\nabla \cdot \vec{u} = 0

b. Navier Stokes equation

\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + \vec{F_B}

c. Energy equation

\rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = - \nabla \cdot (\kappa \nabla T) + 2 \mu {\bf D}:{\bf D}

The compressible and incompressible fluid equations can be summarized as follows.

Governing equation Compressibility 非Compressibility
a.Continuity equation \frac{\partial \rho}{\partial t} + (\vec{u} \cdot \nabla) \rho = - \rho \nabla \cdot \vec{u} $ \nabla \cdot \vec{u} = 0$
b.Navier Stokes equation $\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + (\mu + \lambda) \nabla (\nabla \cdot \vec{u}) + \vec{F_B} $ $\rho \frac{\partial \vec{u}}{\partial t} + \rho (\vec{u} \cdot \nabla)\vec{u} = - \nabla p + \mu \nabla^2 \vec{u} + \vec{F_B} $
c.Energy equation $ \rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = -p \nabla \cdot \vec{u} - \nabla \cdot (\kappa \nabla T) + {\bf T}_{\nu} : {\bf D} + \dot{L}$ $ \rho \frac{\partial e}{\partial t} + \rho (\vec{u} \cdot \nabla) e = - \nabla \cdot (\kappa \nabla T) + 2 \mu {\bf D}:{\bf D}$

2. Incompressible fluid calculation algorithm

Since we are aiming to simulate liquids, this chapter will describe how to solve incompressible fluids. Basically, we will build a discretized equation using an incompressible equation. To avoid complications, we will ignore the energy equations here and mean to solve "explicitly" with the explicit method and "solve implicitly" with the implicit method. Basically, I refer to Compressible Fluid Dynamics.

Typical calculation algorithms are summarized in the table below. The MAC (Maker And Cell) method, Fractional Step method, and SIMPLE method are the basic calculation algorithms for incompressible fluids. It's implemented on various sites and books, so I'll just explain the principle here. At the end of this chapter, we mention the CCUP (CIP-Combined Unified Procedure) method and implement it in the next chapter. The CCUP method is a method for solving the fluid equation by the CIP method used in the previous articles.

algorithm Advection term Diffusion term Pressure term
MAC system Explicit method Explicit method Implicit method
Fractional Step method Explicit method Implicit method Implicit method
---- Implicit method Explicit method Implicit method
SIMPLE method system Implicit method Implicit method Implicit method
CCUP method Explicit method Explicit methodorImplicit method Implicit method

2-1. MAC method

The MAC (Maker And Cell) method is a basic solution method in fluid calculation because it can solve various flows in a stable manner. The calculation procedure is as follows.

  1. Solve the pressure equation by the implicit method to find the pressure $ p ^ {n + 1} $.
  2. Time-advance the velocity equation to find the velocity $ \ vec {u ^ {n + 1}} $.

The pressure equation and the velocity equation are calculated as follows. Multiply the Navier Stokes equation by $ \ nabla $ to discretize and diverge.

\nabla^2 p^{n+1} = - \nabla \cdot \left\{(\vec{u^n} \cdot \nabla )\vec{u^n} \right\} + \frac{\nabla \cdot \vec{u^n}}{\Delta t} + \nabla \cdot \left\{ \nu \nabla^2 \vec{u^n} + \frac{1}{\rho}\vec{F}_B\right\}

It can be obtained. Also, when the advection term and viscosity term of the Navier Stokes equation are explicitly time-discriminated,

    \frac{\vec{u^{n+1}} - \vec{u^n}}{\Delta t} = - (\vec{u} \cdot \nabla)\vec{u^n} - \frac{1}{\rho} \nabla p^n + \nu \nabla^2 \vec{u^n} \\
    \nabla \cdot \vec{u^{n+1}} = 0

It can be obtained. By deliberately solving $ \ nabla \ cdot \ vec {u ^ n} $ without setting it to 0, we are trying to suppress the error during discretization.

There are SMAC (Simplified MAC) method, HSMAC (Highly Simplified MAC) method, etc. as the calculation algorithm of the MAC method system. These were proposed to solve the MAC method at high speed.

2-2. Fractional Step method

In the MAC method, the viscosity term is solved by the explicit method, and it is greatly affected by the time step width limitation (CFL condition). The Fractional Step method described in this section is a method for solving the problem of CFL conditions by implicitly evaluating the viscosity term. Also known as the partial step method. As a feature, in the equation related to velocity used in the MAC method, velocity and pressure were mixed on the right side, but in the Fractional Step method, the velocity and pressure on the right side are completely separated as follows.

\frac{\bar{u} - \vec{u^n}}{\Delta t} = - (\vec{u} \cdot \nabla)\vec{u^n} + \nu \nabla^2 \vec{u^n} \\
    \frac{\vec{u^{n+1}} - \bar{u}}{\Delta t} = - \frac{1}{\rho} \nabla p^{n+1} \\
    \nabla \cdot \vec{u^{n+1}} = 0
\nabla^2 p^{n+1} = \frac{\nabla \cdot \bar{u}}{\Delta t}

As a calculation procedure,

  1. Find the tentative velocity $ \ bar {u} $ from the top equation of the velocity equation.
  2. Use the pressure equation to find the pressure $ p ^ {n + 1} $.
  3. Find the velocity $ \ vec {u ^ {n + 1}} $ using the second equation from the top of the velocity equation.

2-3. SIMPLE method

Abbreviation for Semi-Implicit Method for Pressure-Linked Equation method, which is formulated based on the implicit method. want to know more! Basics of thermo-fluid analysis 64 Chapter 6 Thermo-fluid analysis method: 6.5.4 SIMPLE method explains in detail.

  1. Implicitly find the tentative velocity $ \ bar {u} $.
\bar{u} = \vec{u^n} - \Delta t \left\{(\bar{u} \cdot \nabla)\bar{u} + \nabla p - \nu \nabla^2 \bar{u}\right\}
  1. Update the pressure using the pressure equation.
\nabla^2 \delta p = \frac{\nabla \cdot \bar{u}}{\Delta t}\\ p^{n+1} = p + \delta p
  1. Update the velocity $ \ vec {u} $ with the calculated pressure.
\vec{u^{n+1}} - \bar{u} = - \Delta t \nabla(p^{n+1} - p)

2-4. CCUP method

The CCUP method is a method constructed by Yabe et al., Who devised the CIP method, like when the Burgers equation in previous article was solved by the CIP method. , The governing equation of the fluid is calculated separately for the advection phase and the non-advection phase. Shimizu et al. and Himeno et al. ) Paper is easy to understand. Since the order of solving advection and non-advection terms differs from person to person, we will solve them in the order introduced in the latter paper below. As an algorithm,

  1. In the advection phase, transfer $ density: \ rho $, $ velocity: \ vec {u} $, $ internal energy: e $, $ pressure: p $.
  2. Move to the non-advection phase.
  3. Update velocity, internal energy and pressure at the diffusion stage. (Density is constant)
  4. Solve the pressure equation implicitly and update the pressure.
  5. Find the density, velocity and internal energy using the new pressure.

I will solve in the order. In short, after advection and diffusion of substances, I think that the movement speed is overwhelmingly fast (the same speed as the speed of sound) and the Tsuji 褄 is adjusted.

ccup_method.png

2-5. Spatial difference

Before we get into the implementation, we can implement it without knowing about the method of spatial difference, but I will briefly mention it. This is because pressure vibration called checkerboard instability may occur depending on where the variable is placed in space.

Regular grid

regular_stencil.png

Staggered lattice

staggered_grid.png

Colocate grid

collocated_grid_2.png

In the above method, the basic policy is to use the staggered grid for calculation, but in the implementation of this article, the numerical calculation is performed with the regular grid (~~ I could not understand the motivation when using the staggered grid etc. So ~~).

3. Implementation

Implement by CCUP method. The problem to solve is a cavity flow, which feels like giving a lateral velocity to the upper wall, causing a clockwise wind in the area. The code is posted on Github. Ignore the energy equation. The advection and diffusion terms were mentioned in the previous article [https://qiita.com/KQTS/items/0c4f6c47a4d56881a178), so the new content in the code is the pressure equation part, but with the calculation of the diffusion term. The structure is almost the same. In the pressure calculation, some advection terms are omitted. Also, the boundary conditions are pretty crude

cavity_flow.png

If you express the speed with a red line and the pressure with a blue contour, you can see that a graph like that is output.

2d-cavity.gif

scipy.sparse.dia_matrix

We use scipy.sparse.dia_matrix to create the $ {\ bf A} $ matrix of one-dimensional simultaneous equations $ {\ bf A} \ vec {x} = \ vec {b} $. This matrix $ {\ bf A} $ is called a sparse matrix, and since most of the elements of the matrix are composed of 0, it is specified by scipy.sparse.dia_matrix Diagonal compression storage method (DIA). is efficient in terms of memory and so on. Other storage formats for sparse matrices are described in [Python] Articles that enable high-speed sparse matrix calculations.

An example is shown below.

#Prepare one or more arrays of the same length. Note that you need the same length.
data = np.array([np.arange(10), np.arange(20,30)], np.arange(30,40))  
#Determine the offset of the diagonal matrix according to the index of data. 0:Diagonal matrix,Positive number:Upper side,Negative number:Lower
#Note that if offset is positive, the array will be trimmed from the beginning, and if offset is negative, the array will be trimmed from the back.
offsets = np.array([0, 2, -3])  
a_matrix = scipy.sparse.dia_matrix((data, offsets), shape=(len(data[0]), len(data[0])))  #Set the shape of the matrix with shape.

print(a_matrix)
# <10x10 sparse matrix of type '<class 'numpy.int64'>'
# 	with 25 stored elements (3 diagonals) in DIAgonal format>

print(a_matrix.todense())
# matrix([[ 0,  0, 22,  0,  0,  0,  0,  0,  0,  0],
#         [ 0,  1,  0, 23,  0,  0,  0,  0,  0,  0],
#         [ 0,  0,  2,  0, 24,  0,  0,  0,  0,  0],
#         [30,  0,  0,  3,  0, 25,  0,  0,  0,  0],
#         [ 0, 31,  0,  0,  4,  0, 26,  0,  0,  0],
#         [ 0,  0, 32,  0,  0,  5,  0, 27,  0,  0],
#         [ 0,  0,  0, 33,  0,  0,  6,  0, 28,  0],
#         [ 0,  0,  0,  0, 34,  0,  0,  7,  0, 29],
#         [ 0,  0,  0,  0,  0, 35,  0,  0,  8,  0],
#         [ 0,  0,  0,  0,  0,  0, 36,  0,  0,  9]])

4. Summary

Next, I would like to talk about gas-liquid two-phase (gas and liquid) simulation.

References

Recommended Posts

[Python] Fluid simulation: Incompressible Navier-Stokes equation
[Python] Fluid simulation: Implement the advection equation
[Python] Fluid simulation: Implement the diffusion equation
[Python] Fluid simulation: From linear to non-linear
ambisonics simulation python
First neuron simulation with NEURON + Python
ambisonics simulation (external problem) Python
X-ray diffraction simulation in Python