一维激波管python

一维激波管是一种流体力学问题，可以使用Python进行模拟和可视化。以下是一维激波管的Python实现步骤：

1.导入必要的库和模块

import numpy as np
import matplotlib.pyplot as plt

2.定义初始条件和常数

# 初始条件
rho_L = 1.0  # 左侧密度
u_L = 0.0  # 左侧速度
p_L = 1.0  # 左侧压力
rho_R = 0.125  # 右侧密度
u_R = 0.0  # 右侧速度
p_R = 0.1  # 右侧压力

# 常数
gamma = 1.4  # 比热比

3.定义计算函数

def calc(p):
    """
    计算激波管中的物理量
    """
    if p <= p_R:
        # 右侧区域
        rho = rho_R
        u = u_R
        p = p_R
    elif p > p_R and p <= p_star:
        # 左侧稀疏波区域
        rho = rho_R * ((p / p_R) ** (1 / gamma))
        u = u_R + (p - p_R) * np.sqrt(0.5 * (1 / rho_R + 1 / rho) / p_R)
    elif p > p_star and p <= p_L:
        # 左侧激波区域
        rho = rho_L * ((2 / (gamma + 1) + (gamma - 1) / (gamma + 1) * p / p_L) / (2 * gamma / (gamma + 1) * p / p_L + (gamma - 1) / (gamma + 1))) ** (1 / (gamma - 1))
        u = u_L + 2 * np.sqrt(gamma / (gamma + 1) * p_L / rho_L) * (1 - (p / p_L) ** ((gamma - 1) / (2 * gamma)))
    else:
        # 左侧区域
        rho = rho_L
        u = u_L
        p = p_L
    return rho, u, p

4.定义计算网格和时间步长

# 计算网格和时间步长
dx = 0.01  # 网格间距
x = np.arange(-0.5, 1.5, dx)  # 网格
dt = 0.0002  # 时间步长
t = 0.0  # 初始时间
t_end = 0.2  # 结束时间

5.定义初始状态

# 初始状态
p = np.zeros_like(x)
u = np.zeros_like(x)
rho = np.zeros_like(x)
for i in range(len(x)):
    if x[i] < 0:
        rho[i] = rho_L
        u[i] = u_L
        p[i] = p_L
    else:
        rho[i] = rho_R
        u[i] = u_R
        p[i] = p_R

6.进行计算和可视化

# 进行计算和可视化
while t < t_end:
    # 计算当前时间步长的物理量
    rho_n, u_n, p_n = rho.copy(), u.copy(), p.copy()
    for i in range(1, len(x) - 1):
        rho[i], u[i], p[i] = 0.5 * (rho_n[i - 1] + rho_n[i + 1]) - 0.5 * dt / dx * (rho_n[i + 1] * u_n[i + 1] - rho_n[i - 1] * u_n[i - 1]), \
                             0.5 * (u_n[i - 1] + u_n[i + 1]) - 0.5 * dt / dx * ((rho_n[i + 1] * u_n[i + 1] ** 2 + p_n[i + 1]) - (rho_n[i - 1] * u_n[i - 1] ** 2 + p_n[i - 1])), \
                             0.5 * (p_n[i - 1] + p_n[i + 1]) - 0.5 * dt / dx * (u_n[i + 1] * (rho_n[i + 1] * u_n[i + 1] + p_n[i + 1]) - u_n[i - 1] * (rho_n[i - 1] * u_n[i - 1] + p_n[i - 1]))
    t += dt

# 可视化结果
plt.plot(x, rho, label='Density')
plt.plot(x, u, label='Velocity')
plt.plot(x, p, label='Pressure')
plt.legend()
plt.show()