二维超声速流——普朗特-迈耶稀疏波的流场CFD解（附完整代码）

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

实现了 第 8.3 节 普朗特-迈耶稀疏波流场的数值解  的代码，采用的是 MacCormack 方法，对守恒型方程求解。

注意：

代码最后运行的结果与作者在书中所提供的并不完全相同：

1、当采用作者在书中提供的思路，即在流场下边界壁面处均采用向前差分计算时，结果是在壁面处速度大于稀疏波后速度，稀疏波后速度约710m/s，近壁面速度可达720m/s以上，而作者在书中指出近壁面速度约为706-708m/s。因此在代码中采用了线性外插值的方法。考虑到MacCormack方法采用了两次一阶差分，实际上牺牲了流场上下边界的网格点，因此个人认为最后采用线性插值的方法求解上下边界的网格点更为合理。采用线性插值方法，在近壁面处得到的结果与作者提供的相近。

2、在稀疏波波头处得到的结果光滑，但是在稀疏波波尾处出现震荡（见下图）。修改柯朗数或者人工粘性系数均无法改善。希望有人能提出改进意见。此外，网格数加密了100倍以得到与网格无关的结果。

不足之处，欢迎指正！

运行结果如下：

                                         

   

   

代码如下：

# -*- coding: utf-8 -*-
"""
Created on Sun Aug  9 11:59:45 2020

@author: PANSONG
"""


import numpy as np
import matplotlib.pyplot as plt;plt.close('all')

######################################################
############## 二维超声速流：普朗特-迈耶稀疏波 #########
######################################################

## 几何条件，建模
E = 10 # m, 扩张角位置
theta = 5.352*np.pi/180 # rad，扩张角

## 计算域，流场外形
H = 40 # m, 最大 y 坐标， min_y < 0
L = 65 # m, 最大 x 坐标, min_x = 0

## 指定材料, 物性
gamma = 1.4 # 比热比
R = 1.01e5/1.23/286.1 # 气体常数，p/rho/T

## 边界条件，计算起始条件, 等效初值
Ma1 = 2
p1 = 1.01e5 # N/m2,Pa
T1 = 286.1 # K
rho1 = 1.23 # kg/m3

## 运行参数
C = 0.5 # 柯朗数
Cy = 0.6 # 人工粘性系数
max_Iteration = 10000 # 最大迭代次数

## 离散化， 画网格; 沿 x 方向推进求解，只画 eta 
max_eta = 1  # 计算空间内 eta 最大为 1
num_eta = 401 # eta 方向网格数

eta = np.linspace(0,max_eta,num_eta).reshape(-1,1)

d_eta = max_eta/(num_eta-1)
ksi = [0] #起始计算坐标为 0

y = np.linspace(0,H,num_eta).reshape(-1,1) # 记录网格点 y 坐标

########## 解析解,同时用于校正边界条件 ##########################
def Prandt_Meyer_func(Ma,f2):
    f = np.sqrt((gamma+1)/(gamma-1))*np.arctan(np.sqrt((gamma-1)/(gamma+1)*(Ma**2-1)))-np.arctan(np.sqrt(Ma**2-1))-f2
    return f

## 根据 状态1 和 夹角 得到 状态2 的参数
def Solution_analytic(Ma1,p1,T1,rho1,phi,theta):
    
    f1 = Prandt_Meyer_func(Ma1,0)
    f2 = f1+phi
    
    # 二分法求马赫数
    xm = np.array([1,2.5,5])
    eps = 1e-5
    n = 0
    while abs(xm[0]-xm[2])>= eps:
        fm = Prandt_Meyer_func(xm, f2)
        if fm[1]*fm[2]<0:
            xm = np.array([xm[1],(xm[1]+xm[2])/2,xm[2]])
        else:
            if fm[1]*fm[2]>0:
                xm = np.array([xm[0],(xm[0]+xm[1])/2,xm[1]])
            else:
                xm = np.array([xm[1],xm[1],xm[1]])
        n=n+1
        if n>1000:
            print('Warning ! Max iteration for Ma2')
            break
        
    Ma2 = xm[1]   
    p2 = p1*np.power((1+(gamma-1)/2*Ma1**2)/(1+(gamma-1)/2*Ma2**2),gamma/(gamma-1))
    T2 = T1*(1+(gamma-1)/2*Ma1**2)/(1+(gamma-1)/2*Ma2**2)
    # rho2 = rho1*np.power((1+(gamma-1)/2*Ma1**2)/(1+(gamma-1)/2*Ma2**2),1/(gamma-1))
    rho2 = p2/R/T2
    u2 = Ma2*np.sqrt(gamma*R*T2/(1+np.tan(theta)**2))
    v2 = - u2*np.tan(theta)
    return Ma2,p2,T2,rho2,u2,v2

## 求解析解
Ma2,p2,T2,rho2,u2,v2 = Solution_analytic(Ma1,p1,T1,rho1,theta,theta)

################################################
########## CFD 解 : 守恒型控制方程 #########################
def height(x,eta): #根据输入的位置 x 求对应的 h, 以及 eta 对 x 的偏导数，角度 theta
    # x 是实数，eta 是ndarray
    ## 单位 m, 
    if x <= E:
        h = H
        P_eta_x = np.zeros(shape=eta.shape)
        angle = 0
    else:
        h = H + (x-E)*np.tan(theta)
        P_eta_x = (1-eta)*np.tan(theta)/h
        angle = theta        
    return h,P_eta_x,angle

def calculate_F(rho,u,v,p): # 求F,输入实数或列向量，输出矩阵或行向量
    F1 = rho*u
    F2 = rho*u**2+p
    F3 = rho*u*v
    F4 = gamma/(gamma-1)*p*u+rho*u*(u**2+v**2)/2
    F = np.hstack([F1,F2,F3,F4])
    return F

def calculate_G(F,rho): # 求G,输入F是矩阵，rho是列向量，输出矩阵
    F1 = F[:,0].reshape(-1,1)
    F2 = F[:,1].reshape(-1,1)
    F3 = F[:,2].reshape(-1,1)
    # F4 = F[:,3]
    G1 = rho*F3/F1
    G2 = F3
    G3 = rho*(F3/F1)**2 + F2 - F1**2/rho
    G4 = gamma/(gamma-1)*(F2-F1**2/rho)*F3/F1+rho/2*F3/F1*((F1/rho)**2+(F3/F1)**2)
    G = np.hstack([G1,G2,G3,G4])
    return G

def calculate_original(F): #计算原变量，输入F是矩阵，输出是列向量
    F1 = F[:,0].reshape(-1,1)
    F2 = F[:,1].reshape(-1,1)
    F3 = F[:,2].reshape(-1,1)
    F4 = F[:,3].reshape(-1,1)
    
    A = F3**2/2/F1-F4
    B = gamma/(gamma-1)*F1*F2
    C = -(gamma+1)/2.0/(gamma-1)*F1**3
    rho = (-B + np.sqrt(B**2-4*A*C))/A/2
    
    u = F1/rho
    v = F3/F1
    p = F2-F1*u
    T = p/rho/R
    return rho,u,v,p,T


### 初始化 ######################################
Ma = Ma1*np.ones([num_eta,1])
p = p1*np.ones([num_eta,1])
rho = rho1*np.ones([num_eta,1])
T = T1*np.ones([num_eta,1])
u = Ma*np.sqrt(gamma*R*T) # m/s, x 方向速度
v = np.zeros([num_eta,1]) # m/s, y 方向速度
       
Ma_field = Ma.copy() ## 列数未知，选择在计算中动态增加
p_field = p.copy()
rho_field = rho.copy()
T_field = T.copy()
v_field = v.copy()
u_field = u.copy()

####### MacCormack 方法 ###################################
F = calculate_F(rho, u, v, p)

# 初始化一些要用到的中间变量
P_F_ksi = np.zeros(shape=F.shape)
SF = np.zeros(shape=F.shape)

P_F_ksi_pred = np.zeros(shape=F.shape)
SF_pred = np.zeros(shape=F.shape)

for i in range(max_Iteration):
    
    if i%20==0:
        print('Iteration = ',i)
    
    if ksi[-1] >= L:
        break
    
    ## 计算步长, 坐标变换因子
    hx,P_eta_x,angle = height(ksi[-1],eta) ## x=ksi
    
    dy = d_eta*hx # 物理空间 y 方向网格宽度；网格等宽
    mu = np.arcsin(1/Ma) # 马赫角
    dksi = C*dy/max(max(np.abs(np.tan(theta+mu))),max(np.abs(np.tan(theta+mu))))
    dksi = dksi[0]
    
    y = np.hstack([y,np.linspace(H-hx,H,num_eta).reshape(-1,1)])
    
    ###### 预估步 ################################## 
    
    G = calculate_G(F,rho)
    
    P_F_ksi[:-1,:] = - (P_eta_x[:-1]*(F[1:,:]-F[:-1,:])/d_eta + 1/hx*(G[1:,:]-G[:-1,:])/d_eta)  

    # 人工粘性
    SF[1:-1,:] = Cy*np.abs(p[2:]-2*p[1:-1]+p[:-2])/(p[2:]+2*p[1:-1]+p[:-2])*(F[2:,:]-2*F[1:-1,:]+F[:-2,:])
    SF[0,:] = 2*SF[1,:]-SF[2,:] # 插值壁面的人工粘度
        
    F_pred = F + P_F_ksi*dksi + SF  # 最后一行数据不变，仍是原来的值
    
    # 原变量的预测值
    rho_pred,__,__,p_pred,__ = calculate_original(F_pred)
   
    ###### 校正步 ##############################################
    G_pred = calculate_G(F_pred,rho_pred)
    
    P_F_ksi_pred[1:-1,:] = - (P_eta_x[1:-1]*(F_pred[1:-1,:]-F[:-2,:])/d_eta + 1/hx*(G_pred[1:-1,:]-G_pred[:-2,:])/d_eta)  
       
    P_F_ksi_av = 0.5*(P_F_ksi+P_F_ksi_pred)
    
    #人工粘性
    SF_pred[1:-2,:] = Cy*np.abs(p_pred[2:-1]-2*p_pred[1:-2]+p_pred[:-3])/(p_pred[2:-1]+2*p_pred[1:-2]+p_pred[:-3])*(F_pred[2:-1,:]-2*F_pred[1:-2,:]+F_pred[:-3,:])
    SF_pred[-2,:] = 2*SF_pred[-3,:]-SF_pred[-4,:] #插值上边界倒数第二行的人工粘度
    
    F[1:-1] = F[1:-1] + P_F_ksi_av[1:-1]*dksi + SF_pred[1:-1] ## 校正值,最后一行数据不变；两次差分，损失首尾
    
    
    # 上边界保持不变;最后一行数据不变,无需插值
    # 下边界：先插值，再校正
    F[0,:] = 2*F[1,:]-F[2,:]
    
    # 壁面边界条件 ： 速度与壁面相切; 校正速度方向
    rho,u,v,p,T = calculate_original(F)  #求原变量，都是列向量
    
    phi = angle + np.arctan(v[0]/u[0])
        
    Ma = np.sqrt(v**2+u**2)/np.sqrt(gamma*R*T)

    Ma[0],p[0],T[0],rho[0],u[0],v[0] = Solution_analytic(Ma[0],p[0],T[0],rho[0],phi,angle) # 校正壁面值
    
    F[0,:] = calculate_F(rho[0], u[0], v[0], p[0]) #更新 F
    
    ## 保存所需原变量   
    
    p_field = np.hstack([p_field,p])
    rho_field = np.hstack([rho_field,rho])
    T_field = np.hstack([T_field,T])
    v_field = np.hstack([v_field,v])
    u_field = np.hstack([u_field,u])
    Ma_field = np.hstack([Ma_field,Ma])
    
    ksi.append(ksi[-1]+dksi) #记录每次推进后 ksi 的位置

#%% 计算结果可视化

x = np.tile(ksi,(y.shape[0],1)) ## 扩展成所需 x 的坐标

def mapping(X,Y,Z,name='none'):
    plt.figure()    
    pic = plt.contourf(X,Y,Z,alpha=0.8,cmap='jet')
    plt.colorbar(pic)
    plt.xlabel('x position (m)')
    plt.ylabel('y position (m)')
    plt.title('Evolution of '+name)

mapping(x,y,Ma_field,name='Mach number')
mapping(x,y,rho_field,name='density (kg/m3)')
mapping(x,y,T_field,name='temperature (K)')
mapping(x,y,p_field/1000,name='pressure (kPa)')

length = len(ksi) ## 总 x 数量
x1 = round(length/3)
x2 = round(length/2)
x3 = length - 1
plt.figure()
plt.plot(y[:,x1],u_field[:,x1])
plt.plot(y[:,x2],u_field[:,x2])
plt.plot(y[:,x3],u_field[:,x3])
plt.title('Horizotal velocity in different x positon')
plt.ylabel('Horizontal velocity (m/s)')
plt.xlabel('y position (m)')
plt.legend(['x='+str(round(ksi[x1],1)),'x='+str(round(ksi[x2],1)),'x='+str(round(ksi[x3],1))],loc='upper right')