我需要获得某一时刻的位置速度曲线图像
时间: 2024-03-04 14:49:07 浏览: 29
好的,你可以在代码中添加一些绘图部分,来获得某一时刻的位置速度曲线图像。具体来说,你需要在时间循环中添加一个绘图部分,来绘制某一时刻的速度场和位置速度曲线。以下是修改后的代码:
```python
import numpy as np
from scipy.linalg import solve_banded
import matplotlib.pyplot as plt
# 设置参数
nx = 101 # x方向网格数
ny = 101 # y方向网格数
Lx = 1.0 # 区域长度
Ly = 1.0 # 区域宽度
dx = Lx / (nx - 1) # x方向网格间距
dy = Ly / (ny - 1) # y方向网格间距
dt = 0.001 # 时间步长
T = 0.1 # 计算总时间
nu = 1.0 / 5000 # 动力粘性系数
alpha = dt * nu / dx**2 # 离散化参数
# 初始化网格
x = np.linspace(0, Lx, nx)
y = np.linspace(0, Ly, ny)
X, Y = np.meshgrid(x, y)
u = np.zeros((ny, nx))
v = np.zeros((ny, nx))
p = np.zeros((ny, nx))
# 设定边界条件
u[:, 0] = 0 # 左边界
u[:, -1] = 1 # 右边界
# 设置带状矩阵
a = np.zeros(nx * ny) # 下带
b = np.zeros(nx * ny) # 中带
c = np.zeros(nx * ny) # 上带
d = np.zeros(nx * ny) # 右侧向量
for i in range(ny):
for j in range(nx):
k = i * nx + j
if i == 0 or i == ny - 1 or j == 0 or j == nx - 1:
a[k] = 0
b[k] = 1
c[k] = 0
d[k] = 0 if i == 0 else 1
else:
a[k] = -alpha
b[k] = 1 + 2 * alpha
c[k] = -alpha
d[k] = u[i, j] + dt * (
-0.25 / dx * (u[i, j+1]**2 - u[i, j-1]**2)
-0.25 / dy * (u[i+1, j]*v[i+1, j] - u[i-1, j]*v[i-1, j])
+ nu / dx**2 * (u[i, j+1] - 2*u[i, j] + u[i, j-1])
+ nu / dy**2 * (u[i+1, j] - 2*u[i, j] + u[i-1, j])
)
# 时间循环
t = 0
t_plot = 0.05 # 绘图时刻
u_plot = np.zeros(nx)
while t < T:
# 使用托马斯算法求解带状线性方程组
u_new = np.zeros((ny, nx))
for i in range(ny):
d_new = d[i*nx : (i+1)*nx]
u_new[i] = solve_banded((1, 1), np.vstack((a[i*nx : (i+1)*nx], b[i*nx : (i+1)*nx], c[i*nx : (i+1)*nx])), d_new)
u = u_new
# 更新速度场和压力场
v[1:-1, 1:-1] = v[1:-1, 1:-1] + dt * (
-0.25 / dx * (u[1:-1, 2:] + u[1:-1, 1:-1])**2
+0.25 / dx * (u[1:-1, 1:-1] + u[1:-1, :-2])**2
-0.25 / dy * (u[2:, 1:-1] + u[1:-1, 1:-1])*(v[2:, 1:-1] + v[1:-1, 1:-1])
+0.25 / dy * (u[1:-1, 1:-1] + u[:-2, 1:-1])*(v[1:-1, 1:-1] + v[:-2, 1:-1])
+ nu / dx**2 * (v[1:-1, 2:] - 2*v[1:-1, 1:-1] + v[1:-1, :-2])
+ nu / dy**2 * (v[2:, 1:-1] - 2*v[1:-1, 1:-1] + v[:-2, 1:-1])
)
p[1:-1, 1:-1] = p[1:-1, 1:-1] + dt * (
-0.5 / dx * (u[1:-1, 2:] - u[1:-1, :-2])
-0.5 / dy * (v[2:, 1:-1] - v[:-2, 1:-1])
)
p[0, :] = p[1, :]
p[:, 0] = p[:, 1]
p[-1, :] = p[-2, :]
p[:, -1] = 0
# 更新时间
t = t + dt
# 绘制速度场和位置速度曲线
if t >= t_plot:
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.contourf(X, Y, u, cmap='coolwarm')
plt.colorbar()
plt.title('Velocity Field t={}'.format(t))
plt.xlabel('x')
plt.ylabel('y')
plt.gca().set_aspect('equal')
u_plot = u[int(ny/2), :] # 取中心线的速度
x_plot = x
plt.subplot(122)
plt.plot(x_plot, u_plot)
plt.title('Velocity Profile t={}'.format(t))
plt.xlabel('x')
plt.ylabel('u')
plt.show()
t_plot = t_plot + 0.05 # 每隔0.05秒绘制一次
# 绘制最终的速度场和位置速度曲线
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.contourf(X, Y, u, cmap='coolwarm')
plt.colorbar()
plt.title('Velocity Field t={}'.format(t))
plt.xlabel('x')
plt.ylabel('y')
plt.gca().set_aspect('equal')
plt.subplot(122)
plt.plot(x_plot, u_plot)
plt.title('Velocity Profile t={}'.format(t))
plt.xlabel('x')
plt.ylabel('u')
plt.show()
```
在上面的代码中,我添加了一个 `t_plot` 变量,用来控制绘图的时间。在每次循环中,当时间 `t` 大于等于 `t_plot` 时,就会绘制速度场和位置速度曲线,并将 `t_plot` 增加0.05秒,以便下次绘制。最后,我们再绘制一次最终的速度场和位置速度曲线。
希望这个代码能够帮到你!如果还有其他问题,请继续提问。
相关推荐
![pdf](https://img-home.csdnimg.cn/images/20210720083512.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)