用有限体积方法模拟圆柱绕流的程序

圆柱绕流是一种经典的流体力学问题，可以用有限体积方法来模拟。以下是一个简单的程序示例，使用了基于MacCormack格式的有限体积方法来求解二维不可压缩流体在圆柱周围的流动。

import numpy as np
import matplotlib.pyplot as plt

# 定义常量
L = 1.0  # 长度（正方形域的边长）
N = 100  # 网格数
dx = L / N  # 网格间距
dt = 0.001  # 时间步长
T = 1.0  # 总时间
nu = 0.1  # 动力粘度系数
rho = 1.0  # 密度

# 初始化场变量
u = np.zeros((N+2, N+2))  # x方向速度分量
v = np.zeros((N+2, N+2))  # y方向速度分量
p = np.zeros((N+2, N+2))  # 压力
f = np.zeros((N+2, N+2))  # x方向速度的源项
g = np.zeros((N+2, N+2))  # y方向速度的源项

# 定义边界条件
u[:, 0] = 1.0  # 左边界，入流速度
u[:, -1] = 0.0  # 右边界，自由出流边界
v[0, :] = 0.0  # 上边界，固定壁面边界
v[-1, :] = 0.0  # 下边界，固定壁面边界

# 定义辅助函数
def build_up_b(p, u, v, rho, dt, dx, dy):
    b = np.zeros_like(p)
    b[1:-1, 1:-1] = (rho * (1.0 / dt * ((u[1:-1, 2:] - u[1:-1, :-2]) / (2 * dx) +
                                         (v[2:, 1:-1] - v[:-2, 1:-1]) / (2 * dy)) -
                            ((u[1:-1, 2:] - u[1:-1, :-2]) / (2 * dx))**2 -
                            2 * ((u[2:, 1:-1] - u[:-2, 1:-1]) / (2 * dy) *
                                 (v[1:-1, 2:] - v[1:-1, :-2]) / (2 * dx)) -
                            ((v[2:, 1:-1] - v[:-2, 1:-1]) / (2 * dy))**2))
    return b

def pressure_poisson(p, dx, dy, b):
    pn = np.empty_like(p)
    pn[:] = p[:]
    
    for q in range(100):
        pn[:] = p[:]
        p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, :-2]) * dy**2 +
                           (pn[2:, 1:-1] + pn[:-2, 1:-1]) * dx**2) /
                          (2 * (dx**2 + dy**2)) -
                         dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * b[1:-1, 1:-1])

        p[:, -1] = p[:, -2]  # 右边界，自由出流边界
        p[0, :] = p[1, :]  # 上边界，固定壁面边界
        p[:, 0] = p[:, 1]  # 左边界，入流速度边界
        p[-1, :] = 0  # 下边界，固定壁面边界
    
    return p

# 主循环
for n in range(int(T/dt)):
    un = u.copy()
    vn = v.copy()
    
    # 计算源项
    f[1:-1, 1:-1] = (un[1:-1, 1:-1] -
                     un[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] - un[1:-1, :-2]) -
                     vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[:-2, 1:-1]) -
                     dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, :-2]) +
                     nu * (dt / dx**2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, :-2]) +
                           dt / dy**2 * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[:-2, 1:-1])))
    
    g[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
                     un[1:-1, 1:-1] * dt / dx * (vn[1:-1, 1:-1] - vn[1:-1, :-2]) -
                     vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[:-2, 1:-1]) -
                     dt / (2 * rho * dy) * (p[2:, 1:-1] - p[:-2, 1:-1]) +
                     nu * (dt / dx**2 * (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, :-2]) +
                           dt / dy**2 * (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[:-2, 1:-1])))
    
    # 计算预测速度场
    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, :-2]) -
                     vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[:-2, 1:-1]) -
                     dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, :-2]) +
                     nu * (dt / dx**2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, :-2]) +
                           dt / dy**2 * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[:-2, 1:-1])) +
                     dt * f[1:-1, 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, :-2]) -
                     vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[:-2, 1:-1]) -
                     dt / (2 * rho * dy) * (p[2:, 1:-1] - p[:-2, 1:-1]) +
                     nu * (dt / dx**2 * (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, :-2]) +
                           dt / dy**2 * (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[:-2, 1:-1])) +
                     dt * g[1:-1, 1:-1])
    
    # 处理边界条件
    u[:, 0] = 1.0  # 左边界，入流速度
    u[:, -1] = 0.0  # 右边界，自由出流边界
    v[0, :] = 0.0  # 上边界，固定壁面边界
    v[-1, :] = 0.0  # 下边界，固定壁面边界
    
    # 计算压力场
    b = build_up_b(p, u, v, rho, dt, dx, dy)
    p = pressure_poisson(p, dx, dy, b)
    
    # 打印进度
    if n % 100 == 0:
        print('Iteration: %d' % n)

# 绘制速度场
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.contourf(np.linspace(0, L, N+2), np.linspace(0, L, N+2), np.sqrt(u[1:-1, 1:-1]**2 + v[1:-1, 1:-1]**2))
ax.set_aspect('equal')
plt.show()

在这个程序中，我们使用了MacCormack格式进行时间推进，用中心差分方法离散了空间方程，用SIMPLE算法求解了压力-速度耦合方程。程序的输出是速度场的轮廓图，可以通过调整参数和边界条件来探索不同的圆柱绕流问题。