详细讲解计算流体力学激波,并列出它的方程和代码

计算流体力学（CFD）激波是指在气体中传播的一种压缩波，它是由气体的压缩和加速引起的。激波通常出现在高速流动、喷气发动机、飞行器进出口等应用中，这些应用中需要对气体的流动进行精确的数值模拟。

激波的数学模型可以使用一组非线性偏微分方程来描述，称为欧拉方程组。该方程组包括质量守恒方程、动量守恒方程和能量守恒方程。这些方程可以用以下方式表述：

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v})=0$$

$$\frac{\partial (\rho \mathbf{v})}{\partial t}+\nabla \cdot (\rho \mathbf{v} \mathbf{v}+p\mathbf{I})=0$$

$$\frac{\partial (\rho E)}{\partial t}+\nabla \cdot (\mathbf{v}(\rho E +p))=0$$

其中，$\rho$ 是气体密度，$\mathbf{v}$ 是速度矢量，$p$ 是气体压力，$\mathbf{I}$ 是单位矩阵，$E$ 是气体总能量。

为了求解这些方程，可以使用数值方法，如有限体积法或有限元法。以下是使用 Python 编写的求解欧拉方程组的示例代码：

import numpy as np

# 设置模拟参数
nx = 81
dx = 0.25
dt = 0.0002
gamma = 1.4

# 初始化气体状态
rho_l = 1.0
v_l = 0.0
p_l = 100000.0
rho_r = 0.125
v_r = 0.0
p_r = 10000.0
x_l = 0.5
x_r = 1.0
u = np.zeros((nx, 3))
for i in range(nx):
    if i < nx/2:
        u[i, 0] = rho_l
        u[i, 1] = rho_l * v_l
        u[i, 2] = p_l/(gamma-1) + 0.5*rho_l*v_l**2
    else:
        u[i, 0] = rho_r
        u[i, 1] = rho_r * v_r
        u[i, 2] = p_r/(gamma-1) + 0.5*rho_r*v_r**2

# 定义守恒量到原始量的转换函数
def primitive_variables(u):
    rho = u[:, 0]
    v = u[:, 1] / rho
    p = (gamma-1) * (u[:, 2] - 0.5*rho*v**2)
    return rho, v, p

# 定义原始量到守恒量的转换函数
def conservative_variables(rho, v, p):
    u1 = rho
    u2 = rho * v
    u3 = p/(gamma-1) + 0.5*rho*v**2
    return np.array([u1, u2, u3]).T

# 定义计算通量的函数
def flux(u):
    rho, v, p = primitive_variables(u)
    f1 = u[:, 1]
    f2 = u[:, 1]**2 / u[:, 0] + p
    f3 = u[:, 1]/u[:, 0] * (u[:, 2] + p)
    return np.array([f1, f2, f3]).T

# 定义计算通量梯度的函数
def flux_gradient(u):
    rho, v, p = primitive_variables(u)
    df1_du = np.array([[0, 1, 0]]*u.shape[0])
    df2_du = np.array([[-(gamma-3)/2*v**2, (3-gamma)*v, gamma-1]]*u.shape[0]).T
    df3_du = np.array([[gamma*v*(p/rho+(gamma-1)/2*v**2), -gamma*(p/rho+(gamma-1)*v**2), gamma*(p/rho+(gamma-1)*v**2)]]*u.shape[0]).T
    return np.array([df1_du, df2_du, df3_du])

# 定义计算通量向量的函数
def Riemann_solver(u_l, u_r):
    rho_l, v_l, p_l = primitive_variables(u_l)
    rho_r, v_r, p_r = primitive_variables(u_r)
    c_l = np.sqrt(gamma*p_l/rho_l)
    c_r = np.sqrt(gamma*p_r/rho_r)
    S_l = np.minimum(v_l-c_l, v_r-c_r)
    S_r = np.maximum(v_l+c_l, v_r+c_r)
    if S_l >= 0:
        f_l = flux(u_l)
        return f_l
    elif S_r <= 0:
        f_r = flux(u_r)
        return f_r
    else:
        f_l = flux(u_l)
        f_r = flux(u_r)
        return (S_r*f_l - S_l*f_r + S_l*S_r*(u_r - u_l)) / (S_r - S_l)

# 定义计算时间步长的函数
def compute_dt(u, dx):
    rho, v, p = primitive_variables(u)
    c = np.sqrt(gamma*p/rho)
    u_abs = np.abs(v) + c
    return 0.9*dx / np.max(u_abs)

# 执行主程序
t = 0.0
while t < 0.01:
    dt = compute_dt(u, dx)
    if t+dt > 0.01:
        dt = 0.01 - t
    u_old = u.copy()
    for i in range(1, nx-1):
        u_l = u_old[i-1, :]
        u_r = u_old[i, :]
        f = Riemann_solver(u_l, u_r)
        u[i, :] = u_old[i, :] - dt/dx * (f[i,:] - f[i-1,:])
    t += dt

# 输出结果
rho, v, p = primitive_variables(u)
x = np.linspace(0, nx*dx, nx)
import matplotlib.pyplot as plt
plt.plot(x, rho, 'r-', label='Density')
plt.plot(x, v, 'b-', label='Velocity')
plt.plot(x, p, 'g-', label='Pressure')
plt.legend()
plt.show()

在上面的代码中，我们使用有限体积法来求解欧拉方程组，并使用 Riemann 求解器来计算通量向量。该代码模拟了一个激波在气体中传播的过程，并输出了气体密度、速度和压力随时间和位置的变化。