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

最新推荐文章于 2023-10-24 12:14:44 发布

Patrick_PS

最新推荐文章于 2023-10-24 12:14:44 发布

阅读量3.9k

点赞数 12

分类专栏： CFD学习笔记文章标签： python cfd

本文链接：https://blog.csdn.net/weixin_45205765/article/details/107888990

版权

CFD学习笔记专栏收录该内容

6 篇文章 22 订阅

订阅专栏

入门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 root

def 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<1
    else:
        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]) #对于给定的 T0

p_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])/dx

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

def calculate_F(U1,U2,U3):
    F1 = U2
    F2 = U2**2/U1+(gamma-1)/gamma*(U3-gamma/2*U2**2/U1)
    F3 = gamma*U2*U3/U1-gamma*(gamma-1)/2*U2**3/U1**2
    F = np.vstack([F1,F2,F3])
    return F

for 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*T
    S = 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*S
        
    U_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 + J
    
    P_U_t_av = 0.5*(P_U_t_pred + P_U_t)
    
    #人工粘性
    rho_pred = U_pred[0,:]/A
    T_pred = (gamma-1)*(U_pred[2,:]/U_pred[0,:]-gamma/2*(U_pred[1,:]/U_pred[0,:])**2)
    p_pred = rho_pred*T_pred
    S_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_pred
    
    U = 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，插值入流边界是为了记录history
    U[2,0] = U[0,0]*(1/(gamma-1)+gamma/2*(U[1,0]/U[0,0])**2) #T[0]=1
 
    rho = U[0,:]/A
    rho_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 root

def 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<1
    else:
        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]) #对于给定的 T0

p_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]) ## CFD
    plt.plot(x[1],y[1]) ## analytic
    plt.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()

Patrick_PS

关注

12
点赞
踩
49

收藏

觉得还不错? 一键收藏
0
评论
喷管流动的守恒型CFD解法及激波捕捉（附完整代码）

入门CFD，主要参考书目《计算流体力学基础及其应用》（John D.Anderson著，吴颂平等译）实现了第 7.6节激波捕捉的代码，采用的是MacCormack 方法，对守恒型方程求解；关于非守恒型方程，可见亚声速-超声速等熵喷管拟一维流动的CFD解法（附完整代码）。代码中增加了求解析解的功能。对求解析解过程的理解和思考如下：1、激波的存在可看做间断点，将激波前后的流动分别看做是两个等熵流动。如下图所示（图片截取自参考书，下标带0的压力表示总压，亦即初始压力）。激波前的等熵流动...
复制链接

扫一扫