别再被湍流模型搞晕了！用Python从零实现一个超简单的DNS求解器（附完整代码）

用Python从零实现极简DNS求解器：让Navier-Stokes方程看得见摸得着

当第一次听说"直接数值模拟"(DNS)时，我盯着那组复杂的Navier-Stokes方程看了整整一个下午——那些偏微分符号像天书一样令人望而生畏。直到有一天，我决定用Python把这些方程变成可以运行的代码，才发现原来理解流体力学可以如此直观。本文将带你用不到200行Python代码，构建一个能真实运行的二维方腔流DNS求解器，让你亲手"触摸"流体运动的每一个细节。

1. 准备工作：理解极简版NS方程

在开始编码前，我们需要对Navier-Stokes方程做适当简化。考虑二维不可压缩流动，控制方程可表示为：

# 连续性方程 ∂u/∂x + ∂v/∂y = 0 # x方向动量方程 ∂u/∂t + u∂u/∂x + v∂u/∂y = -1/ρ ∂p/∂x + ν(∂²u/∂x² + ∂²u/∂y²) # y方向动量方程 ∂v/∂t + u∂v/∂x + v∂v/∂y = -1/ρ ∂p/∂y + ν(∂²v/∂x² + ∂²v/∂y²)

其中：

• u,v分别是x,y方向速度分量
• p是压力
• ρ是流体密度(设为1)
• ν是运动粘度系数

提示：方腔流问题指在一个矩形区域内，顶部壁面以恒定速度移动，其余三面为静止壁面，是验证CFD方法的经典案例。

2. 数值方法：有限差分法实现

2.1 网格生成与变量排列

我们采用交错网格(staggered grid)布置变量，避免压力振荡问题：

import numpy as np # 参数设置 nx, ny = 41, 41 # 网格数 dx = 1.0 / (nx-1) dy = 1.0 / (ny-1) # 初始化变量 (在交错网格位置) u = np.zeros((ny, nx)) # x速度 (位于网格右面) v = np.zeros((ny, nx)) # y速度 (位于网格上面) p = np.zeros((ny, nx)) # 压力 (位于网格中心)

2.2 压力泊松方程求解

压力通过求解泊松方程获得，我们使用松弛迭代法：

def solve_pressure(p, u, v, dx, dy, rho=1.0, nit=50): pn = np.empty_like(p) for _ in range(nit): pn = p.copy() p[1:-1,1:-1] = ((dy**2*(pn[1:-1,2:]+pn[1:-1,0:-2]) + dx**2*(pn[2:,1:-1]+pn[0:-2,1:-1])) - rho*dx**2*dy**2*((u[1:-1,2:]-u[1:-1,0:-2])/(2*dx) + (v[2:,1:-1]-v[0:-2,1:-1])/(2*dy))) / (2*(dx**2+dy**2)) # 边界条件 p[0,:] = p[1,:] # dp/dy = 0 at y=0 p[-1,:] = 0 # p=0 at y=1 p[:,0] = p[:,1] # dp/dx = 0 at x=0 p[:,-1] = p[:,-2] # dp/dx = 0 at x=1 return p

2.3 时间推进算法

采用分步法(time-splitting)进行时间积分：

def time_step(u, v, p, dx, dy, dt, nu=0.1, rho=1.0): # 临时速度场 un = u.copy() vn = v.copy() # 对流项计算 u[1:-1,1:-1] = (un[1:-1,1:-1] - dt/dx*un[1:-1,1:-1]*(un[1:-1,1:-1]-un[1:-1,0:-2]) - dt/dy*vn[1:-1,1:-1]*(un[1:-1,1:-1]-un[0:-2,1:-1]) + 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] - dt/dx*un[1:-1,1:-1]*(vn[1:-1,1:-1]-vn[1:-1,0:-2]) - dt/dy*vn[1:-1,1:-1]*(vn[1:-1,1:-1]-vn[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]))) # 压力修正 p = solve_pressure(p, u, v, dx, dy, rho) # 速度修正 u[1:-1,1:-1] -= dt/dx * (p[1:-1,2:]-p[1:-1,0:-2])/2 v[1:-1,1:-1] -= dt/dy * (p[2:,1:-1]-p[0:-2,1:-1])/2 # 边界条件 u[0,:] = 0 # 下壁面 u[-1,:] = 1 # 上壁面移动 u[:,0] = 0 # 左壁面 u[:,-1] = 0 # 右壁面 v[0,:] = 0 # 下壁面 v[-1,:] = 0 # 上壁面 v[:,0] = 0 # 左壁面 v[:,-1] = 0 # 右壁面 return u, v, p

3. 完整求解流程实现

现在我们将所有部分组合起来，形成完整的求解器：

import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 参数设置 nx, ny = 41, 41 nt = 500 # 时间步数 dt = 0.001 nu = 0.1 # 粘度系数 rho = 1.0 dx = 1.0/(nx-1) dy = 1.0/(ny-1) # 初始化 u = np.zeros((ny,nx)) v = np.zeros((ny,nx)) p = np.zeros((ny,nx)) # 主循环 for n in range(nt): u, v, p = time_step(u, v, p, dx, dy, dt, nu, rho) # 每50步可视化 if n % 50 == 0: fig = plt.figure(figsize=(11,7)) plt.contourf(np.linspace(0,1,nx), np.linspace(0,1,ny), p, alpha=0.5) plt.colorbar() plt.quiver(np.linspace(0,1,nx), np.linspace(0,1,ny), u, v) plt.xlabel('X') plt.ylabel('Y') plt.title(f'Time Step {n}') plt.show()

4. 结果分析与可视化

运行上述代码后，我们可以观察到方腔内逐渐形成的涡旋结构。为了更专业地分析结果，添加流线可视化：

def plot_streamlines(u, v, p): x = np.linspace(0, 1, nx) y = np.linspace(0, 1, ny) X, Y = np.meshgrid(x, y) plt.figure(figsize=(11,7)) plt.streamplot(X, Y, u, v, density=2, color='k', linewidth=1) plt.contourf(X, Y, p, alpha=0.5) plt.colorbar() plt.xlabel('X') plt.ylabel('Y') plt.title('Streamlines and Pressure Contour') plt.show() plot_streamlines(u, v, p)

典型输出应显示：

1. 顶部移动壁面驱动的顺时针主涡
2. 左下角和右下角的次级涡
3. 压力在角落区域达到最大值

5. 性能优化与扩展建议

虽然我们的求解器已经可以工作，但在实际应用中还需要考虑以下优化：

计算效率提升技巧：

• 使用Numba加速数值计算
• 实现多核并行计算
• 采用更高效的线性代数求解器

物理模型扩展方向：

• 添加温度场模拟自然对流
• 实现三维版本求解器
• 引入自适应网格细化(AMR)技术

# Numba加速示例 from numba import jit @jit(nopython=True) def time_step_accelerated(u, v, p, dx, dy, dt, nu, rho): # 实现与之前相同的逻辑，但会被编译为机器码 ...

第一次运行这个求解器时，我惊讶地发现原来那些抽象的数学符号真的能转化为肉眼可见的流体运动。当第一个涡旋在屏幕上出现时，那种"啊哈时刻"的兴奋感至今难忘。建议你在修改参数(如粘度系数、网格分辨率)后重新运行，观察流动模式的变化——这才是理解CFD最有趣的方式。