基于二维欧拉方程、空间格式选用MUSCL格式，时间格式选择一阶前向差分，写一段关于Rayleigh-Tarlor不稳定性流动的python代码

由于Rayleigh-Taylor不稳定性流动是一个三维问题，本代码只实现了其二维模拟，并采用了有限体积法求解。

代码如下：

import numpy as np
import matplotlib.pyplot as plt

# 定义常数
gamma = 5/3        # 气体比热容
rho1 = 1           # 上层密度
rho2 = 2           # 下层密度
P1 = 1             # 上层压力
P2 = 2             # 下层压力
g = 1              # 重力加速度
Lx = 1             # x方向长度
Ly = 1             # y方向长度
Nx = 100           # x方向网格数
Ny = 100           # y方向网格数
dx = Lx / Nx       # x方向网格间距
dy = Ly / Ny       # y方向网格间距
cfl = 0.5          # CFL数
t_final = 1        # 终止时间

# 定义初始条件
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y)
u1 = np.ones((Ny, Nx)) * rho1
u2 = np.zeros((Ny, Nx))
u2[Y > 0.5] = rho2
u = np.stack([u1, u2], axis=-1)
p = np.ones((Ny, Nx)) * P1
p[Y > 0.5] = P2
e = p / ((gamma - 1) * u[..., 0]) + 0.5 * u[..., 1]**2 / u[..., 0]
U = np.concatenate([u, e[..., None]], axis=-1)

# 定义MUSCL格式
def muscl(r):
    r_l = np.roll(r, 1, axis=0)
    r_r = np.roll(r, -1, axis=0)
    r_ll = np.roll(r, 2, axis=0)
    r_rr = np.roll(r, -2, axis=0)
    Dp = (r - r_l) / 2
    Dm = (r_r - r) / 2
    D = np.where(Dp * Dm > 0, np.where(np.abs(Dp) < np.abs(Dm), Dp, Dm), 0)
    r_pl = r + D / 2
    r_mi = r - D / 2
    r_p = (3 * r - 4 * r_l + r_ll) / 2
    r_m = (3 * r_r - 4 * r - r_rr) / 2
    return r_mi, r_pl, r_m, r_p

# 定义数值通量
def numerical_flux(r_l, r_r):
    rho_l, u_l, p_l = r_l[..., 0], r_l[..., 1] / r_l[..., 0], (gamma - 1) * (r_l[..., 2] - 0.5 * r_l[..., 1]**2 / r_l[..., 0])
    rho_r, u_r, p_r = r_r[..., 0], r_r[..., 1] / r_r[..., 0], (gamma - 1) * (r_r[..., 2] - 0.5 * r_r[..., 1]**2 / r_r[..., 0])
    a_l, a_r = np.sqrt(gamma * p_l / rho_l), np.sqrt(gamma * p_r / rho_r)
    H_l, H_r = (r_l[..., 2] + p_l) / rho_l + 0.5 * u_l**2, (r_r[..., 2] + p_r) / rho_r + 0.5 * u_r**2
    S_l, S_r = np.minimum(u_l - a_l, u_r - a_r), np.maximum(u_l + a_l, u_r + a_r)
    F_l, F_r = np.stack([r_l[..., 1], r_l[..., 1] * u_l + p_l, r_l[..., 2] * u_l + p_l * u_l], axis=-1), np.stack([r_r[..., 1], r_r[..., 1] * u_r + p_r, r_r[..., 2] * u_r + p_r * u_r], axis=-1)
    F_star = (S_r * F_l - S_l * F_r + S_l * S_r * (r_r - r_l)) / (S_r - S_l)
    return np.where(S_l > 0, F_l, np.where(S_r < 0, F_r, np.where(S_l * S_r > 0, F_star, 0)))

# 定义二维欧拉方程
def euler_equation(U):
    r = U[..., :2]
    F = np.zeros_like(U)
    F[..., :2] = U[..., 1:]
    F[..., 1] += (gamma - 1) * (U[..., 2] - 0.5 * U[..., 1]**2 / U[..., 0])
    F[..., 2] += gamma * U[..., 2] * (U[..., 1] / U[..., 0]) + (gamma - 1) * U[..., 1] * (U[..., 2] - 0.5 * U[..., 1]**2 / U[..., 0]) / U[..., 0]
    r_mi_x, r_pl_x, r_m_x, r_p_x = muscl(r)
    r_mi_y, r_pl_y, r_m_y, r_p_y = muscl(np.moveaxis(r, [0, 1], [1, 0]))
    F_x = numerical_flux(r_m_x, r_p_x)
    F_y = np.moveaxis(numerical_flux(np.moveaxis(r_m_y, [0, 1], [1, 0]), np.moveaxis(r_p_y, [0, 1], [1, 0])), [1, 0, 2], [0, 1, 2])
    G = F_x - np.roll(F_x, -1, axis=0) + F_y - np.roll(F_y, -1, axis=1)
    return G / (dx * dy)

# 定义前向差分方法
def forward_diff(U, dt):
    return U + dt * euler_equation(U)

# 定义计算CFL数的函数
def calc_cfl(U):
    u = U[..., 1] / U[..., 0]
    a = np.sqrt(gamma * (gamma - 1) * (U[..., 2] / U[..., 0] - 0.5 * u**2))
    return np.max(np.abs(u) + a) * dt / dx

# 主程序
t = 0
while t < t_final:
    dt = cfl * dx / calc_cfl(U)
    U = forward_diff(U, dt)
    t += dt

# 绘图
fig, ax = plt.subplots()
ax.contourf(X, Y, U[..., 0], cmap='jet')
plt.show()

值得注意的是，本代码并没有实现边界条件，因此模拟结果并不准确。实际应用中需要根据具体情况选择合适的边界条件。