python解N-S方程

一、引言
看着手中花花绿绿的CFD软件，再回首先辈留下的旷世秘籍——N-S方程，如同“白月光”一样，照射在我们的床头。我们究竟是离梦想越来越近，还是与梦想渐行渐远了呢？[声针之家微信公众号]
NS方程的重要性不言而喻，他支撑起了流体力学和计算流体力学的学术体系，被称为流体工作者的白月光，笔者偶尔发现了网络上计算NS方程的Python代码，分享给大家，仅供大家参考。
二、代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

# NACA 0012 airfoil function
def naca0012(x):
    return 0.6 * (0.2969 * np.sqrt(x) - 0.1260 * x - 0.3516 * x**2 + 0.2843 * x**3 - 0.1015 * x**4)

# Grid parameters
nx, ny = 100, 100
dx, dy = 1.0 / (nx - 1), 1.0 / (ny - 1)
nt = 500
nit = 50  # Define the number of iterations for the Poisson solver
dt = 0.0001  # Reduced time step to avoid overflow
rho = 1.0
nu = 0.1

# Generate grid
x = np.linspace(0, 1, nx)
y = np.linspace(-0.5, 0.5, ny)
X, Y = np.meshgrid(x, y)

# Initialize fields
u = np.ones((ny, nx))
v = np.zeros((ny, nx))
p = np.zeros((ny, nx))
b = np.zeros((ny, nx))

# Initialize airfoil boundary
airfoil = (Y < naca0012(X)) & (Y > -naca0012(X))

# Function to build the source term b
def build_up_b(b, rho, dt, u, v, dx, dy):
    b[1:-1, 1:-1] = (rho * (1 / dt * ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) +
                                      (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
                           ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx)) ** 2 -
                           2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
                                (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx)) -
                           ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) ** 2))

# Function to solve the Poisson equation for pressure
def pressure_poisson(p, dx, dy, b):
    pn = np.empty_like(p)
    for q in range(nit):
        pn = p.copy()
        p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +
                          (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /
                         (2 * (dx**2 + dy**2)) -
                         dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * b[1:-1, 1:-1])

        # Boundary conditions for pressure
        p[:, -1] = p[:, -2]
        p[0, :] = p[1, :]
        p[:, 0] = p[:, 1]
        p[-1, :] = 0
    return p

# Function to solve the Navier-Stokes equations
def navier_stokes(nt, u, v, dt, dx, dy, p, rho, nu):
    un = np.empty_like(u)
    vn = np.empty_like(v)
    b = np.zeros((ny, nx))

    for n in range(nt):
        un = u.copy()
        vn = v.copy()

        build_up_b(b, rho, dt, un, vn, dx, dy)
        p = pressure_poisson(p, dx, dy, b)

        u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
                         un[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
                         vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
                         dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, 0:-2]) +
                         nu * (dt / dx**2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
                               dt / dy**2 * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])))

        v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
                         un[1:-1, 1:-1] * dt / dx * (vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
                         vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
                         dt / (2 * rho * dy) * (p[2:, 1:-1] - p[0:-2, 1:-1]) +
                         nu * (dt / dx**2 * (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
                               dt / dy**2 * (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))

        # Apply boundary conditions
        u[0, :] = 1  # Inlet
        u[-1, :] = 1  # Outlet
        u[:, 0] = 1  # Inlet
        u[:, -1] = 1  # Outlet
        v[0, :] = 0
        v[-1, :] = 0
        v[:, 0] = 0
        v[:, -1] = 0

        # Apply airfoil boundary conditions
        u[airfoil] = 0
        v[airfoil] = 0

    return u, v, p

# Solve the equations
u, v, p = navier_stokes(nt, u, v, dt, dx, dy, p, rho, nu)

# Plot the results
plt.figure(figsize=(11, 7), dpi=100)
plt.contourf(X, Y, p, alpha=0.5, cmap=plt.cm.viridis)
plt.colorbar()
plt.contour(X, Y, p, cmap=plt.cm.viridis)
plt.streamplot(X, Y, u, v, color='k')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Flow around a NACA 0012 airfoil')
plt.show()

以上是一个机翼绕流的代码，从结果看相比较各形各色的CFD计算软件出的图片差距还是比较大。对于复杂流场自己编程计算结果肯定会出现各种各样的问题，对于现成的CFD软件内部加入了很多平滑机制，使计算结果即使快要发散的情况下尽量能够继续求解下去，毕竟经常崩溃的软件用户体验会很差。因此，很多CFD软件使其数据稳定的方式也是核心技术。