用有限体积法模拟二维圆柱绕流的程序
时间: 2023-12-12 19:03:30 浏览: 359
以下是一个基于有限体积法的二维圆柱绕流模拟程序的示例代码。该程序使用了Python语言和NumPy库。
```python
import numpy as np
# Constants
nx = 101 # Number of grid points in x direction
ny = 101 # Number of grid points in y direction
nt = 1000 # Number of time steps
dt = 0.001 # Time step size
dx = 2 / (nx - 1) # Grid spacing in x direction
dy = 2 / (ny - 1) # Grid spacing in y direction
rho = 1 # Density
nu = 0.1 # Kinematic viscosity
sigma = 0.2 # Relaxation parameter
# Initialization
u = np.zeros((ny, nx)) # Velocity in x direction
v = np.zeros((ny, nx)) # Velocity in y direction
p = np.zeros((ny, nx)) # Pressure
b = np.zeros((ny, nx)) # Source term
# Define function to calculate pressure
def pressure_poisson(p, b, dx, dy):
pn = np.empty_like(p)
pn[:] = p[:]
for q in range(nt):
pn[:] = p[:]
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)
p[1:-1, 1:-1] = (1 - sigma) * p[1:-1, 1:-1] + \
sigma / (2 * (dx ** 2 + dy ** 2)) * \
((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy ** 2 +
(pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx ** 2 +
b[1:-1, 1:-1] * dx ** 2 * dy ** 2)
p[:, -1] = p[:, -2] # Boundary condition: p = 0 at x = 2
p[0, :] = p[1, :] # Boundary condition: dp/dy = 0 at y = 0
p[-1, :] = p[-2, :] # Boundary condition: dp/dy = 0 at y = 2
p[50, 50] = 0 # Pressure at center of cylinder is 0
return p
# Define function to calculate velocity
def cavity_flow(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[:]
vn[:] = v[:]
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)
p = pressure_poisson(p, b, dx, dy)
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])))
u[0, :] = 0 # Boundary condition: u = 0 at y = 0
u[:, 0] = 0 # Boundary condition: u = 0 at x = 0
u[:, -1] = 0 # Boundary condition: u = 0 at x = 2
u[-1, :] = 1 # Boundary condition: u = 1 at y = 2
v[0, :] = 0 # Boundary condition: v = 0 at y = 0
v[:, 0] = 0 # Boundary condition: v = 0 at x = 0
v[:, -1] = 0 # Boundary condition: v = 0 at x = 2
v[-1, :] = 0 # Boundary condition: v = 0 at y = 2
return u, v, p
# Run simulation
u, v, p = cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)
# Plot results
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(np.linspace(0, 2, nx), np.linspace(0, 2, ny))
ax.plot_surface(X, Y, p[:], cmap='viridis')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$');
```
该程序使用了矩阵运算,因此可以相对快速地模拟二维圆柱绕流。您可以根据需要进行修改和调整,以满足特定的模拟需求。
阅读全文