二维超声速流——普朗特-迈耶稀疏波的流场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] #起始计算坐标为 0y = 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))-f2return 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-5n = 0while 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+1if n>1000:print('Warning ! Max iteration for Ma2')breakMa2 = 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/T2u2 = 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 = HP_eta_x = np.zeros(shape=eta.shape)angle = 0else:h = H + (x-E)*np.tan(theta)P_eta_x = (1-eta)*np.tan(theta)/hangle = theta        return h,P_eta_x,angledef calculate_F(rho,u,v,p): # 求F,输入实数或列向量，输出矩阵或行向量F1 = rho*uF2 = rho*u**2+pF3 = rho*u*vF4 = gamma/(gamma-1)*p*u+rho*u*(u**2+v**2)/2F = np.hstack([F1,F2,F3,F4])return Fdef 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/F1G2 = F3G3 = rho*(F3/F1)**2 + F2 - F1**2/rhoG4 = 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 Gdef 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-F4B = gamma/(gamma-1)*F1*F2C = -(gamma+1)/2.0/(gamma-1)*F1**3rho = (-B + np.sqrt(B**2-4*A*C))/A/2u = F1/rhov = F3/F1p = F2-F1*uT = p/rho/Rreturn 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=ksidy = 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')

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

普朗特迈耶稀疏波_埃里克·迈耶的访谈
普朗特迈耶稀疏波 I've always wanted to interview Eric Meyer. His early CSS books are a big reason this blog ...
基于MATLAB的二维与三维插值拟合运算（附完整代码）
· 一. 一维插值 interp1函数在上个博客中(如下链接)已经更新了,此处再补充两个相关例题. 基于MATLAB的数据插值运算:Lagrange与Hermite算法(附完整代码)_唠嗑!的博客-C ...
快速集成二维码扫描，使用最新版本的zxing（2017.11.10抽取zxing代码）
github 地址: github.com/maning0303/- ZXingCode 快速集成二维码扫描,使用最新版本的zxing代码提取(2017.11.10) 功能: 1:生成二维码(带Log ...
OpenCV差分二值化的实时场景文本检测的实例(附完整代码)
OpenCV差分二值化的实时场景文本检测的实例 OpenCV差分二值化的实时场景文本检测的实例 OpenCV差分二值化的实时场景文本检测的实例 OpenCV差分二值化的实时场景文本检测的实例(附完整代 ...
Three.js实例详解___旋转的精灵女孩(附完整代码和资源)（二）
Three.js实例详解___旋转的精灵女孩(附完整代码和资源)(二) 本篇目录: 五.实例中所使用的代码语法详细解释 (1).构建一个三维空间场景 (2).选择一个透视投影相机作为观察点 (a).创 ...
java 稀疏数组和二维数组转换，并保存稀疏数组到文件后可以读取
稀疏数组和二维数组转换稀疏数组:当一个数组中大部分元素为0,或者为同一个值的数组时,可以使用稀疏数组来保存该数组稀疏数组的处理方法: 记录数组一共有多少行,有多少个不同的值把具有不同值得元素的行 ...
python二维向量运算模拟_Python数学基础之向量定义与向量运算（附代码）
患难与困苦是磨练人格的最高学府.--苏格拉底(公元前470年-公元前399年) Adversity and pinch are the highest institution of higher le ...
asp生成带参数的二维码并合成推广海报图片，asp合并合成推广海报图片asp代码
最近做的一个项目中,客户要求用asp生成二维码,然后合并到一张背景图片上,合并生成一张推广海报来,可把我愁坏了,经过一个晚上的努力,成功了,下面把这个:asp生成带参数的二维码并合成推广海报,asp合 ...
asp生成带参数的二维码并合成推广海报图片，asp合并合成推广海报图片asp代码...
最近做的一个项目中,客户要求用asp生成二维码,然后合并到一张背景图片上,合并生成一张推广海报来,可把我愁坏了,经过一个晚上的努力,成功了,下面把这个:asp生成带参数的二维码并合成推广海报,asp合 ...

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

注意：

运行结果如下：

代码如下：

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

最新文章

热门文章