喷管流动的守恒型CFD解法及激波捕捉（附完整代码）

入门CFD，主要参考书目《计算流体力学基础及其应用》（John D.Anderson 著，吴颂平等译）

实现了 第 7.6 节激波捕捉 的代码，采用的是 MacCormack 方法，对守恒型方程求解；关于非守恒型方程，可见亚声速-超声速等熵喷管拟一维流动的CFD解法（附完整代码）。

代码中增加了求解析解的功能。对求解析解过程的理解和思考如下：

1、激波的存在可看做间断点，将激波前后的流动分别看做是两个等熵流动。如下图所示（图片截取自参考书，下标带0的压力表示总压，亦即初始压力）。激波前的等熵流动初始温度为 T0，初始压力为 p0，无初速度；激波后的等熵流动初始温度也为 T0，初始压力为 0.6882p0，无初速度。（系数0.6882的计算请往下看，至于为什么可以认为初温相同，欢迎补充！）

2、关联激波前后等熵的流动的要素是质量守恒定律。质量流量由下式计算：

式中，A* 代表声速喉道面积。在给定喷管入口条件和喷管形状时，如果流动能够达到声速，则 A* 等于实际喉道面积；如果流动均为亚声速，则 A* 小于实际喉道面积，对应的质量流量也小。由 1 可知质量流量与 $p_{0}A^*$ 成正比。由质量守恒定律可以得到 $p_{01}A_1^* = p_{02}A_2^*$ 。

3、结合喷管等熵流动的公式：

可得：

$\frac{p_e}{p_{02}}*\frac{A_e}{A_2^*} = \frac{p_e}{p_{01}}*\frac{A_e}{A_1^*} = f(Ma_e)$

在激波存在时，流动能达到音速， $A_1^*$ 为实际喉道面积， $p_e$ / $p_{01}$ 和 $A_e$ / $A_1^*$ 均可由已知条件确定，据此可以解出 $Ma_e$ 。解出 $Ma_e$ 后，可以计算出 $p_e$ / $p_{02}$ 。最后，根据下述关系，结合已知条件，即可算出 $p_{02}=0.6882p_{01}$ 。

$\frac{p_{02}}{p_{01}}=\frac{p_{02}}{p_e} * \frac{p_e}{p_{01}}$

$p_{02}$ / $p_{01}$ 是波前马赫数的函数，据此可以解出波前马赫数，从而解出该马赫数对应的截面积，得到激波所在位置。

回顾质量守恒，可以得出 $A_2^*$ / $A_1^*$ 的值。最后指出，在同一个问题，统一以激波前的参数为基准确定无量纲物理量的具体数值，因此在计算激波后流动时，需要根据 $A_2^*$ / $A_1^*$ 和 $p_{02}$ / $p_{01}$ 进行换算。

不足之处，欢迎指正！

求解结果：

求解析解的代码：

#################### 解析解 ########################################################
from scipy.optimize import rootdef func(Ma,arg):return (1.0/1.2*(1.0+0.2*Ma**2.0))**3-arg*Ma #马赫数是面积（位置）的函数## 激波前
Ma_analytic_1 = []
x_shock = 2.1 ## 激波位置
for x0 in x[np.where(x<=x_shock)]:A0 = section(x0)r = root(func, x0=[0,10],args=A0)['x'] #有两个，分别对应喉口两侧if x0<1.5:Ma = min(r) # 当 x<1.5 时，Ma<1else:Ma = max(r)Ma_analytic_1.append(Ma)Ma_analytic_1 = np.array(Ma_analytic_1)
## 激波后
p02_p01 = 0.6882 #激波前后总压比
A1star_a2star = p02_p01 #根据质量守恒，激波前后的声速喉口面积比
Ma_analytic_2 = []
for x0 in x[np.where(x>x_shock)]:A0 = section(x0)*A1star_a2star ##  激波后喷管无量纲面积, 分母为 A2*r = root(func, x0=[0,10],args=A0)['x'] #有两个，分别对应喉口两侧Ma = min(r) # 亚声速Ma_analytic_2.append(Ma)Ma_analytic_2 = np.array(Ma_analytic_2)Ma_analytic = np.hstack([Ma_analytic_1,Ma_analytic_2])## 激波后的流动等效为 初始温度也是 T0, 初始压力为 p02 = p0*0.6882 的等熵流动
rho_analytic_1 = np.power(1+(gamma-1)/2*Ma_analytic_1**2,-1/(gamma-1))
rho_analytic_2 = np.power(1+(gamma-1)/2*Ma_analytic_2**2,-1/(gamma-1))
rho_analytic = np.hstack([rho_analytic_1, rho_analytic_2*0.6882]) # p = rho*T, rho 跟 p 情况相同T_analytic_1 = np.power(1+(gamma-1)/2*Ma_analytic_1**2,-1)
T_analytic_2 = np.power(1+(gamma-1)/2*Ma_analytic_2**2,-1)
T_analytic = np.hstack([T_analytic_1, T_analytic_2]) #对于给定的 T0p_analytic_1 = np.power(1+(gamma-1)/2*Ma_analytic_1**2,-gamma/(gamma-1))
p_analytic_2 = np.power(1+(gamma-1)/2*Ma_analytic_2**2,-gamma/(gamma-1))
p_analytic = np.hstack([p_analytic_1, p_analytic_2*0.6882]) # p02/p0 = 0.6882, 统一换成 p0 的无量纲量mass_analytic = rho_analytic*A*Ma_analytic*np.sqrt(T_analytic)####################################################################################

完整代码：

# -*- coding: utf-8 -*-
"""
Created on Wed Aug  5 19:11:15 2020@author: PANSONG
"""#### Laval 喷管 拟一维计算 ###############import numpy as np
import matplotlib.pyplot as plt;plt.close('all')### 守恒型控制方程，激波捕捉
#%% 建模
#喷管截面形状, 几何建模
def section(x):return 1+2.2*(x-1.5)**2# 初始条件
def initial_condition(x):rho1 = 1 + 0*x[np.where(x<=0.5)]    rho2 = 1.0-0.366*(x[np.where((0.5<x)&(x<=1.5))]-0.5)rho3 = 0.634-0.702*(x[np.where((1.5<x)&(x<=2.1))]-1.5)rho4 = 0.5892-0.10228*(x[np.where((2.1<x)&(x<=3.0))]-2.1)T1 = 1.0 + 0*x[np.where(x<=0.5)]T2 = 1.0-0.167*(x[np.where((0.5<x)&(x<=1.5))]-0.5)T3 = 0.833-0.4908*(x[np.where((1.5<x)&(x<=2.1))]-1.5)T4 = 0.93968-0.0622*(x[np.where((2.1<x)&(x<=3.0))]-2.1)rho = np.hstack([rho1,rho2,rho3,rho4])T = np.hstack([T1,T2,T3,T4])V = 0.59/(rho*section(x))return rho,T,V# 边界条件
pe = 0.6784  #出流压力# 理想气体比热比， 指定材料
gamma = 1.4## 离散化， 画网格
max_x = 3
dx = 0.05
num_point = int(max_x/dx + 1)
x = np.linspace(0,3,num_point)mid = np.where(x==1.5)[0] ## 中间点max_Iteration = 1400#%% 动态绘图
plt.ion() #开启interactive mode 成功的关键函数
plt.figure()
plt.title('Dimensionless Parameters Evolution')
plt.xlabel('Iteration')iterations = []#%%  CFD 计算
#### 初始化，无量纲量
A = section(x)## 这里给两个初值，一个 0，一个线性，避免迭代过程中 残差判断时 重复判断 i
array0 = np.ones(shape=x.shape)*0.01 rho,T,V = initial_condition(x)
rho_history = np.vstack([array0,rho])
T_history = np.vstack([array0,T])
V_history = np.vstack([array0,V])## 守恒变量初值计算，这里选择不保存守恒变量
U1 = rho*A
U2 = rho*A*V
U3 = rho*(T/(gamma-1)+gamma/2*V**2)*A
U = np.vstack([U1,U2,U3])# 辅助变量
J = np.zeros(shape=U[:,2:].shape) # J 维度比 U 和 F 少 2
zero = np.zeros(shape=U.shape[0]) # zero 用于填充偏导数矩阵和人工粘性矩阵，3行1列C = 0.5 # 柯朗数
Cx = 0.2 # 人工粘性系数
####### MacCormack 方法 ###################################
# 向前，向后差分计算，输入 N 个数据，输出 N-2 个偏导数
def forward_partial(y,dx):return (y[2:]-y[1:-1])/dxdef backward_partial(y,dx):return (y[1:-1]-y[:-2])/dxdef calculate_F(U1,U2,U3):F1 = U2F2 = U2**2/U1+(gamma-1)/gamma*(U3-gamma/2*U2**2/U1)F3 = gamma*U2*U3/U1-gamma*(gamma-1)/2*U2**3/U1**2F = np.vstack([F1,F2,F3])return Ffor i in range(max_Iteration):# 绘制 喉口处 无量纲温度 和 无量纲密度iterations.append(i)if i%20 == 0:plt.plot(iterations,rho_history[1:,mid],'r-',linewidth=1.0)plt.plot(iterations,T_history[1:,mid],'b-',linewidth=1.0)plt.legend(['Density','Temperature'])plt.pause(0.01)##########################################   a = np.sqrt(T) #无量纲当地声速        dt = min(C*dx/(V[1:-1]+a[1:-1])) # 计算最大允许时间推进步长  rho_error = max(np.abs((rho_history[i]-rho_history[i-1])/dt))V_error = max(np.abs((V_history[i]-V_history[i-1])/dt))T_error = max(np.abs((T_history[i]-T_history[i-1])/dt))error = max(rho_error,V_error,T_error) ## 取最大残差if error < 1e-6:break###### 预估步 ##################################F = calculate_F(U[0,:],U[1,:],U[2,:])# J 维度比 U 和 F 少 2；J[1,:] = (gamma-1)/gamma/A[1:-1]*(U[2,1:-1]-gamma/2*U[1,1:-1]**2/U[0,1:-1])*forward_partial(A, dx)# J[1,:] = 1/gamma*rho[1:-1]*T[1:-1]*forward_partial(A, dx)P_U_t = -(F[:,2:]-F[:,1:-1])/dx + J  # 向量计算# 人工粘性p = rho*TS = np.abs(p[2:]-2*p[1:-1]+p[:-2])/(p[2:]+2*p[1:-1]+p[:-2])*(U[:,2:]-2*U[:,1:-1]+U[:,:-2])S = Cx*SU_pred = U + np.column_stack([zero,P_U_t,zero])*dt + np.column_stack([zero,S,zero])# 完善入口出预估值U_pred[1,0] = 2*U_pred[1,1]-U_pred[1,2] # U2, 预估步只插值入流边界，出流边界未用在校正步上U_pred[2,0] = U_pred[0,0]*(1/(gamma-1)+gamma/2*(U_pred[1,0]/U_pred[0,0])**2) # U3###### 校正步 ##############################################F_pred = calculate_F(U_pred[0,:],U_pred[1,:],U_pred[2,:])J[1,:] = (gamma-1)/gamma/A[1:-1]*(U_pred[2,1:-1]-gamma/2*U_pred[1,1:-1]**2/U_pred[0,1:-1])*backward_partial(A, dx)# J_pred[1,:] = 1/gamma*rho_pred[1:-1]*T_pred[1:-1]*backward_partial(A, dx)P_U_t_pred = -(F_pred[:,1:-1]-F_pred[:,:-2])/dx + JP_U_t_av = 0.5*(P_U_t_pred + P_U_t)#人工粘性rho_pred = U_pred[0,:]/AT_pred = (gamma-1)*(U_pred[2,:]/U_pred[0,:]-gamma/2*(U_pred[1,:]/U_pred[0,:])**2)p_pred = rho_pred*T_predS_pred = np.abs(p_pred[2:]-2*p_pred[1:-1]+p_pred[:-2])/(p_pred[2:]+2*p_pred[1:-1]+p_pred[:-2])*(U_pred[:,2:]-2*U_pred[:,1:-1]+U_pred[:,:-2])S_pred = Cx*S_predU = U + np.column_stack([zero,P_U_t_av,zero])*dt + np.column_stack([zero,S_pred,zero])# 完善出口校正值U[:,-1] = 2*U[:,-2]-U[:,-3] # 校正步插值出流边界，用在下一次迭代的预估步VN = U[1,-1]/U[0,-1]U[2,-1] = pe*A[-1]/(gamma-1)+gamma/2*U[1,-1]*VN #保证出口处的压力边界条件## 计算原变量   U[1,0] = 2*U[1,1]-U[1,2] # U2，插值入流边界是为了记录historyU[2,0] = U[0,0]*(1/(gamma-1)+gamma/2*(U[1,0]/U[0,0])**2) #T[0]=1rho = U[0,:]/Arho_history = np.vstack([rho_history,rho])V = U[1,:]/U[0,:]V_history = np.vstack([V_history,V])T = (gamma-1)*(U[2,:]/U[0,:]-gamma/2*V**2)T_history = np.vstack([T_history,T])plt.ioff() #关闭动态绘图#%% 计算结果可视化
#################### 解析解 ########################################################
from scipy.optimize import rootdef func(Ma,arg):return (1.0/1.2*(1.0+0.2*Ma**2.0))**3-arg*Ma #马赫数是面积（位置）的函数## 激波前
Ma_analytic_1 = []
x_shock = 2.1 ## 激波位置
for x0 in x[np.where(x<=x_shock)]:A0 = section(x0)r = root(func, x0=[0,10],args=A0)['x'] #有两个，分别对应喉口两侧if x0<1.5:Ma = min(r) # 当 x<1.5 时，Ma<1else:Ma = max(r)Ma_analytic_1.append(Ma)Ma_analytic_1 = np.array(Ma_analytic_1)
## 激波后
p02_p01 = 0.6882 #激波前后总压比
A1star_a2star = p02_p01 #根据质量守恒，激波前后的声速喉口面积比
Ma_analytic_2 = []
for x0 in x[np.where(x>x_shock)]:A0 = section(x0)*A1star_a2star ##  激波后喷管无量纲面积, 分母为 A2*r = root(func, x0=[0,10],args=A0)['x'] #有两个，分别对应喉口两侧Ma = min(r) # 亚声速Ma_analytic_2.append(Ma)Ma_analytic_2 = np.array(Ma_analytic_2)Ma_analytic = np.hstack([Ma_analytic_1,Ma_analytic_2])## 激波后的流动等效为 初始温度也是 T0, 初始压力为 p02 = p0*0.6882 的等熵流动
rho_analytic_1 = np.power(1+(gamma-1)/2*Ma_analytic_1**2,-1/(gamma-1))
rho_analytic_2 = np.power(1+(gamma-1)/2*Ma_analytic_2**2,-1/(gamma-1))
rho_analytic = np.hstack([rho_analytic_1, rho_analytic_2*0.6882]) # p = rho*T, rho 跟 p 情况相同T_analytic_1 = np.power(1+(gamma-1)/2*Ma_analytic_1**2,-1)
T_analytic_2 = np.power(1+(gamma-1)/2*Ma_analytic_2**2,-1)
T_analytic = np.hstack([T_analytic_1, T_analytic_2]) #对于给定的 T0p_analytic_1 = np.power(1+(gamma-1)/2*Ma_analytic_1**2,-gamma/(gamma-1))
p_analytic_2 = np.power(1+(gamma-1)/2*Ma_analytic_2**2,-gamma/(gamma-1))
p_analytic = np.hstack([p_analytic_1, p_analytic_2*0.6882]) # p02/p0 = 0.6882, 统一换成 p0 的无量纲量mass_analytic = rho_analytic*A*Ma_analytic*np.sqrt(T_analytic)####################################################################################
def plot_results(x,y,x0,y0,title='value'):plt.figure()plt.plot(x[0],y[0]) ## CFDplt.plot(x[1],y[1]) ## analyticplt.scatter(x0,y0,color='red')plt.xlabel('Dimensionless distance')plt.title('Steady '+title+' in different position')plt.legend(['CFD solution','Analytic solution','throat'])# plot_results([x,x],[rho_history[-1,:],rho_analytic],1.5,rho_history[-1,mid],title='density')
# plot_results([x,x],[T_history[-1,:],T_analytic],1.5,T_history[-1,mid],title='temperature')p_history = rho_history*T_history
plot_results([x,x],[p_history[-1,:],p_analytic],1.5,p_history[-1,mid],title='pression')Ma_history = V_history/np.sqrt(T_history)
plot_results([x,x],[Ma_history[-1,:],Ma_analytic],1.5,Ma_history[-1,mid],title='Mach number')mass_history = rho_history*A*V_history
plot_results([x,x],[mass_history[-1,:],mass_analytic],1.5,mass_history[-1,mid],title='mass flow')
plt.ylim([0.4,0.7])plt.show()