高超声速钝头体外流代码示例

以下是一份 Python 代码示例，用于计算高超声速钝头体外流：

import numpy as np
import matplotlib.pyplot as plt

# 定义常数
p1 = 1.225  # 空气密度
p2 = 1.4*p1  # 燃气密度
R = 287.1  # 气体常数
T1 = 300  # 空气温度
a1 = np.sqrt(1.4*R*T1)  # 空气声速
u1 = 1500  # 入口速度
M1 = u1/a1  # 入口马赫数
theta = 20  # 钝头半角
L = 0.2  # 钝头长度
L1 = L/np.sin(np.deg2rad(theta))  # 圆锥长
L2 = L*np.tan(np.deg2rad(theta))  # 钝头长
A1 = np.pi/4*(0.5**2)  # 入口面积
A2 = np.pi/4*(0.5**2)  # 出口面积

# 计算出口速度和压强
M2 = np.sqrt((M1**2 + 2/(1.4-1)*(1 - (p2/p1)*(A1/A2)*(1 + (1.4/2)*(M1**2*np.sin(np.deg2rad(theta))**2 - 1))))  # 出口马赫数
u2 = M2*a1*np.sin(np.deg2rad(theta))  # 出口速度
T2 = T1*(1 + 0.5*(1.4-1)*M1**2*np.sin(np.deg2rad(theta))**2)  # 出口温度
p2 = p1*(1 + 0.5*(1.4-1)*M1**2*np.sin(np.deg2rad(theta))**2)**(1.4/(1.4-1))  # 出口压强

# 计算流场参数
x = np.linspace(0, L1+L2, 1000)
p = np.zeros_like(x)
T = np.zeros_like(x)
rho = np.zeros_like(x)
u = np.zeros_like(x)
M = np.zeros_like(x)
p[0] = p1
T[0] = T1
rho[0] = p1/(R*T1)
u[0] = u1
M[0] = M1
for i in range(1, len(x)):
    if x[i] <= L1:
        p[i] = p[i-1]*(1 + 0.5*(1.4-1)*M[i-1]**2)**(-1.4/(1.4-1))
        T[i] = T[i-1]*(1 + 0.5*(1.4-1)*M[i-1]**2)**(-1)
        rho[i] = p[i]/(R*T[i])
        u[i] = M[i-1]*a1*np.sin(np.deg2rad(theta))
        M[i] = np.sqrt((M[i-1]**2 + 2/(1.4-1)*(1 - (p[i]/p[i-1])*(A1/A2)*(1 + (1.4/2)*(M[i-1]**2*np.sin(np.deg2rad(theta))**2 - 1)))))
    else:
        p[i] = p2
        T[i] = T2
        rho[i] = p2/(R*T2)
        u[i] = u2
        M[i] = M2

# 绘制图像
plt.plot(x, p/p1, label='Pressure ratio')
plt.plot(x, T/T1, label='Temperature ratio')
plt.plot(x, rho/rho[0], label='Density ratio')
plt.plot(x, u/a1, label='Velocity ratio')
plt.plot(x, M, label='Mach number')
plt.xlabel('Position along flow field (m)')
plt.legend()
plt.show()

这段代码计算了一个高超声速钝头体的出口速度、压强和流场参数，然后用 Matplotlib 将结果绘制成图像。请注意，这只是一份示例代码，具体计算需要根据具体问题进行调整和修改。