亚声速-超声速等熵喷管拟一维流动的CFD解法（附完整代码）

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

实现了 第 7.3 节 亚声速-超声速等熵喷管流动的CFD解法 的代码，采用的是 MacCormack 方法。代码中增加了动态显示无量纲温度和无量纲密度的功能，参考的是 python中plot实现即时数据动态显示方法，最终结果如下：图中红点代表喉口所在位置。

不足之处，欢迎指正 ！

   

        

                                

完整代码如下：

# -*- coding: utf-8 -*-
"""
Created on Wed Aug  5 19:11:15 2020

@author: PANSONG
"""

#### Laval 喷管 拟一维计算 ###############

import numpy as np
import matplotlib.pyplot as plt

#%% 参数，采用无量纲量
#喷管截面形状
def section(x):
    return 1+2.2*(x-1.5)**2

# 初始条件
def initial_condition(x):
    rho = 1-0.314*x
    T = 1-0.2314*x
    V = (0.1+1.09*x)*np.sqrt(T)
    return rho,T,V

# 理想气体比热比
gamma = 1.4

## 离散化
max_x = 3 
dx = 0.1
num_point = int(max_x/dx + 1)
x = np.linspace(0,3,num_point)

max_Iteration = 10000

#%% 动态绘图
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])

####### MacCormack 方法 ###################################
# 向前，向后差分计算，输入 N 个数据，输出 N-2 个偏导数
def forward_partial(y,dx):
    return (y[2:]-y[1:-1])/dx

def backward_partial(y,dx):
    return (y[1:-1]-y[:-2])/dx

def new_rho_T(original,partial,dt):
    new = original + np.hstack([0,partial,0])*dt # P_rho_t 首尾补 0 
    new[-1] = 2*new[-2]-new[-3]  # 出口处又外插值确定，入口处保持不变
    return new

def new_V(V,P_V_t,dt):
    new_V = V + np.hstack([0,P_V_t,0])*dt
    new_V[0] = 2*new_V[1] - new_V[2] #入口的速度可变，外插值确定
    new_V[-1] = 2*new_V[-2]-new_V[-3]
    return new_V

C = 1 # 柯朗数,小于等于 1;从精度考虑，尽可能接近 1
# 关于柯朗数，通过对线性双曲型偏微分方程的稳定性分析，要求 C<=1，对应数值解的依赖区域包含精确解依赖区域
# 从解的精度考虑，应使数值解的依赖区域尽可能等于精确解依赖区域，即 C 尽可能接近 1
# 对于非线性偏微分方程，起指导性作用，参见参考书目 P221

for i in range(max_Iteration):
    
    # 绘制 喉口处 无量纲温度 和 无量纲密度
    iterations.append(i)
    plt.plot(iterations,rho_history[1:,15],'r-',linewidth=1.0)
    plt.plot(iterations,T_history[1:,15],'b-',linewidth=1.0)
    plt.legend(['Density','Temperature'])
    plt.pause(1e-3)
    
    ##########################################   
    rho_error = max(np.abs((rho_history[i]-rho_history[i-1])/rho_history[i-1]))
    V_error = max(np.abs((V_history[i]-V_history[i-1])/V_history[i-1]))
    T_error = max(np.abs((T_history[i]-T_history[i-1])/T_history[i-1]))
    error = max(rho_error,V_error,T_error) ## 取最大相对误差
    if error < 1e-4:
        break
    
    a = np.sqrt(T) #无量纲当地声速        
    dt = min(C*dx/(V+a)) # 计算最大允许时间推进步长   
    ###### 预估步 ########
    P_rho_t = -V[1:-1]*forward_partial(rho, dx) - rho[1:-1]*forward_partial(V, dx) - rho[1:-1]*V[1:-1]*forward_partial(np.log(A), dx)
    
    rho_pred = new_rho_T(rho,P_rho_t,dt)
    
    P_V_t = -V[1:-1]*forward_partial(V, dx) - (1.0/gamma)*(forward_partial(T, dx)+(T[1:-1]/rho[1:-1])*forward_partial(rho, dx))
    
    V_pred = new_V(V,P_V_t,dt)
    
    P_T_t = -V[1:-1]*forward_partial(T, dx) - (gamma-1.0)*T[1:-1]*(forward_partial(V, dx)+V[1:-1]*forward_partial(np.log(A), dx))
    
    T_pred = new_rho_T(T,P_T_t,dt)
    
    ###### 校正步 ########
    P_rho_t_pred = ( -V_pred[1:-1]*backward_partial(rho_pred, dx) 
                   -rho_pred[1:-1]*backward_partial(V_pred, dx) 
                   -rho_pred[1:-1]*V_pred[1:-1]*backward_partial(np.log(A), dx) )
    
    P_rho_t_av = 0.5*( P_rho_t_pred + P_rho_t )
    
    P_V_t_pred = ( -V_pred[1:-1]*backward_partial(V_pred, dx)
                   -1.0/gamma*backward_partial(T_pred, dx)
                   -1.0/gamma*T_pred[1:-1]/rho_pred[1:-1]*backward_partial(rho_pred, dx) )
    
    P_V_t_av = 0.5*( P_V_t_pred + P_V_t )
    
    P_T_t_pred = ( -V_pred[1:-1]*backward_partial(T_pred, dx)
             -(gamma-1.0)*T_pred[1:-1]*backward_partial(V_pred, dx)
             -(gamma-1.0)*T_pred[1:-1]*V_pred[1:-1]*backward_partial(np.log(A), dx) )
    
    P_T_t_av = 0.5*( P_T_t_pred + P_T_t )
    
    rho = new_rho_T(rho,P_rho_t_av,dt)
    rho_history = np.vstack([rho_history,rho])
    
    V = new_V(V,P_V_t_av,dt)
    V_history = np.vstack([V_history,V])
    
    T = new_rho_T(T,P_T_t_av,dt)
    T_history = np.vstack([T_history,T])
    
plt.ioff() #关闭动态绘图

#%% 计算结果可视化
def plot_results(x,y,x0,y0,title='value'):
    
    plt.figure()
    plt.plot(x,y)
    plt.scatter(x0,y0,color='red')
    plt.xlabel('Dimensionless distance')
    plt.title('Steady '+title+' in different position')

plot_results(x,rho_history[-1,:],1.5,rho_history[-1,15],title='density')
plot_results(x,T_history[-1,:],1.5,T_history[-1,15],title='temperature')

p = rho_history[-1,:]*T_history[-1,:]
plot_results(x,p,1.5,p[15],title='pression')

Ma = V_history[-1,:]/np.sqrt(T_history[-1,:])
plot_results(x,Ma,1.5,Ma[15],title='Mach number')

plt.show()