计算流体力学（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 编写的求解欧拉方程组的示例代码：
```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 求解器来计算通量向量。该代码模拟了一个激波在气体中传播的过程，并输出了气体密度、速度和压力随时间和位置的变化。
3、能量方程的特殊形式
在这里回答一下我上次的疑问，就是总压在可压缩流动中的含义，在这里我们可以看到，所有滞止参数的表达式都是从能量方程推导得到，之前不可压缩的总压定义在这里已经不在适用，总压是根据等熵关系式得到，并不能像伯努利方程中的总压一样代表着机械能，这里的理解可能不太严谨，后续可能会有更好的理解。
286
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
