|
入门CFD，主要参考书目《计算流体力学基础及其应用》（John D.Anderson 著，吴颂平等 译）
实现了 第 10 章 流过平板的超声速流动 的代码。利用有限差分法求解二维 Navier-Stokes 方程，采用MacCormack 显示方法。代码总体不难，按照作者提供的思路编写即可。最后得到的结果中，压力一项与作者给出的有出入，在平板后缘靠近壁面的地方有震荡，猜测与边界条件有关：书中给出的壁面边界条件是对压力进行插值计算，而代码中选择了水平速度，垂直速度，密度和温度作为独立变量进行计算，在壁面处对密度进行插值。未对此猜测进行进一步验证。其他变量，如温度，速度，马赫数等得到的结果与作者提供的基本相同。
不足之处，欢迎指正 ！
结果如下：
 
  
 
 
完整代码如下：
# -*- coding: utf-8 -*-
"""
Created on Sun Aug 16 14:49:14 2020
@author: PANSONG
"""
import numpy as np
import matplotlib.pyplot as plt;plt.close('all')
######################################################
############## 流过平板的超声速流动 ###################
######################################################
## 几何条件，建模，## 计算域，流场外形
L = 0.00001 # m, 平板长度
## 指定材料, 物性
gamma = 1.4 # 比热比
R = 287 # J/kg/K 气体常数
Pr = 0.71 # 普朗特常数
cv = R/(gamma-1) # 定容比热
cp = cv + R # 定压比热
mu0 = 1.789e-5 # Pa·s 参考动力粘度
T0 = 288.16 # K， 参考温度
#### 边界条件
# 来流
Ma_inf = 4
a_inf = 340.28 # m/s
p_inf = 101325.0 # Pa
T_inf = 288.16 # K
v_inf = 0 # 速度垂直分量
u_inf = Ma_inf*a_inf # 速度水平分量
ReL = p_inf/R/T_inf*u_inf*L/mu0 # 来流雷诺数
delta = 5*L/np.sqrt(ReL) # 边界层厚度
H = 5*delta # 计算域高度
rho_inf = p_inf/R/T_inf
# 恒温壁面
Tw = 1*T_inf
## 运行参数
K = 0.6 # 柯朗数
max_Iteration = 8000 # 最大迭代次数
## 离散化， 画网格;
xnum = 71 # 计算空间内 eta 最大为 1
ynum = 71 # eta 方向网格点数为
x = np.linspace(0,L,xnum)
y = np.linspace(0,H,ynum)
dx = L/(xnum-1)
dy = H/(ynum-1)
###########################################################
#### CFD 求解 : 守恒型控制方程, 二维 Navier-Stokes 方程 #####
###########################################################
############# 初始化 ######################################
u = u_inf*np.ones([xnum,ynum]) # 选择 u,v,rho,T 作为独立变量
v = v_inf*np.ones([xnum,ynum])
rho = rho_inf*np.ones([xnum,ynum])
T = T_inf*np.ones([xnum,ynum])
u[:,0] = 0
T[1:,0] = Tw
u_history = np.zeros([max_Iteration,xnum,ynum])
v_history = np.zeros([max_Iteration,xnum,ynum])
rho_history = np.zeros([max_Iteration,xnum,ynum])
T_history = np.zeros([max_Iteration,xnum,ynum])
# 初始化一些要用到的中间变量
U = np.zeros([4,xnum,ynum])
E = np.zeros([4,xnum,ynum])
F = np.zeros([4,xnum,ynum])
U_pred = 0.001*np.ones([4,xnum,ynum]) # 避免出现除以 0
E_pred = np.zeros([4,xnum,ynum])
F_pred = np.zeros([4,xnum,ynum])
def tau_xx(mu,u,v,case=0): # 计算粘性应力
tau_xx = np.zeros(shape=u.shape)
P_v_y = np.zeros(shape=u.shape)
P_u_x = np.zeros(shape=u.shape)
P_v_y[:,1:-1] = (v[:,2:]-2*v[:,1:-1]+v[:,:-2])/dy/2.0 # tau_xx 只出现在 E里，对 y 中心差分
P_v_y[:,0] = (v[:,1]-v[:,0])/dy # 物面只能单侧差分
P_v_y[:,-1] = (v[:,-1]-v[:,-2])/dy
if case == 0:
P_u_x[:-1,:] = (u[1:,:]-u[:-1,:])/dx # 向前差分
P_u_x[-1,:] = (u[-1,:]-u[-2,:])/dx # 只能单侧差分，相当于零阶外插
tau_xx = -2.0/3.0*mu*(P_u_x+P_v_y)+2*mu*P_u_x
elif case == 1:
P_u_x[1:,:] = (u[1,:]-u[:-1,:])/dx # 向后差分
P_u_x[0,:] = (u[1,:]-u[0,:])/dx
tau_xx = -2.0/3.0*mu*(P_u_x+P_v_y)+2*mu*P_u_x
else:
tau_xx = None
print('Wrong assignment of case in calculate tau_xx')
return tau_xx
def tau_yy(mu,u,v,case=0): # 计算粘性应力
tau_yy = np.zeros(shape=u.shape)
P_v_y = np.zeros(shape=u.shape)
P_u_x = np.zeros(shape=u.shape)
P_u_x[1:-1,:] = (u[2:,:]-2*u[1:-1,:]+u[:-2,:])/dx/2.0 # tau_yy 只出现在 F里，对 x 中心差分
P_u_x[0,:] = (u[1,:]-u[0,:])/dx
P_u_x[-1,:] = (u[-1,:]-u[-2,:])/dx
if case == 0:
P_v_y[:,:-1] = (v[:,1:]-v[:,:-1])/dy # 向前差分
P_v_y[:,-1] = (v[:,-1]-v[:,-2])/dy
tau_yy = -2.0/3.0*mu*(P_v_y+P_u_x)+2*mu*P_v_y
elif case == 1:
P_v_y[:,1:] = (v[:,1:]-v[:,:-1])/dy # 向后差分
P_v_y[:,0] = (v[:,1]-v[:,0])/dy
tau_yy = -2.0/3.0*mu*(P_v_y+P_u_x)+2*mu*P_v_y
else:
tau_yy = None
print('Wrong assignment of case in calculate tau_yy')
return tau_yy
def tau_xy_E(mu,u,v,case=0): # 计算粘性应力
tau_xy = np.zeros(shape=u.shape)
P_u_y = np.zeros(shape=u.shape)
P_v_x = np.zeros(shape=u.shape)
P_u_y[:,1:-1] = (u[:,2:]-2*u[:,1:-1]+u[:,:-2])/dy/2.0 # 对 y 中心差分
P_u_y[:,0] = (u[:,1]-u[:,0])/dy
P_u_y[:,-1] = (u[:,-1]-u[:,-2])/dy
if case == 0:
P_v_x[:-1,:] = (v[1:,:]-v[:-1,:])/dx # 向前差分
P_v_x[-1,:] = (v[-1,:]-v[-2,:])/dx
tau_xy = mu*(P_v_x+P_u_y)
elif case == 1:
P_v_x[1:,:] = (v[1:,:]-v[:-1,:])/dx # 向后差分
P_v_x[0,:] = (v[1,:]-v[0,:])/dx
tau_xy = mu*(P_v_x+P_u_y)
else:
tau_xy = None
print('Wrong assignment of case in calculate tau_xx_E')
return tau_xy
def tau_xy_F(mu,u,v,case=0): # 计算粘性应力
tau_xy = np.zeros(shape=u.shape)
P_u_y = np.zeros(shape=u.shape)
P_v_x = np.zeros(shape=u.shape)
P_v_x[1:-1,:] = (v[2:,:]-2*v[1:-1,:]+v[:-2,:])/dx/2.0 # 对 x 中心差分
P_v_x[0,:] = (v[1,:]-v[0,:])/dx
P_v_x[-1,:] = (v[-1,:]-v[-2,:])/dx
if case == 0:
P_u_y[:,:-1] = (u[:,1:]-u[:,:-1])/dy # 向前差分
P_u_y[:,-1] = (u[:,-1]-u[:,-2])/dy
tau_xy = mu*(P_u_y+P_v_x)
elif case == 1:
P_u_y[:,1:] = (u[:,1:]-u[:,:-1])/dy # 向后差分
P_u_y[:,0] = (u[:,1]-u[:,0])/dy # 相当于零阶外插
tau_xy = mu*(P_u_y+P_v_x)
else:
tau_xy = None
print('Wrong assignment of case in calculate tau_xx_F')
return tau_xy
def qx(T,k,case=0): # 计算粘性应力
qx = np.zeros([xnum,ynum])
if case == 0:
qx[:-1,:] = -k[:-1,:]*(T[1:,:]-T[:-1,:])/dx # 向前差分
qx[-1,:] = qx[-2,:]
elif case == 1:
qx[1:,:] = -k[1:,:]*(T[1:,:]-T[:-1,:])/dx # 向后差分
qx[0,:] = qx[1,:]
else:
qx = None
print('Wrong assignment of case in calculate qx')
return qx
def qy(T,k,case=0): # 计算粘性应力
qy = np.zeros([xnum,ynum])
if case == 0:
qy[:,:-1] = -k[:,:-1]*(T[:,1:]-T[:,:-1])/dy # 向前差分
qy[:,-1] = qy[:,-2]
elif case == 1:
qy[:,1:] = -k[:,1:]*(T[:,1:]-T[:,:-1])/dy # 向后差分
qy[:,0] = qy[:,1]
else:
qy = None
print('Wrong assignment of case in calculate qy')
return qy
def calculate_U(rho,u,v,e): #
U1 = rho
U2 = rho*u
U3 = rho*v
U5 = rho*(e+(u**2+v**2)/2)
U = np.array([U1,U2,U3,U5])
return U
def calculate_original(U): #计算原变量
rho = U[0,:,:]
u = U[1,:,:]/U[0,:,:]
v = U[2,:,:]/U[0,:,:]
e = U[3,:,:]/U[0,:,:]-(u**2+v**2)/2
T = e/cv
return rho,u,v,T
def calculate_E(rho,u,v,e,p,tau_xx,tau_xy,qx): #
E1 = rho*u
E2 = rho*u**2 + p - tau_xx
E3 = rho*u*v - tau_xy
E5 = (rho*(e+(u**2+v**2)/2) + p)*u - u*tau_xx - v*tau_xy + qx
E = np.array([E1,E2,E3,E5])
return E
def calculate_F(rho,u,v,e,p,tau_yy,tau_xy,qy): #
F1 = rho*v
F2 = rho*u*v - tau_xy
F3 = rho*v**2 + p - tau_yy
F5 = (rho*(e+(u**2+v**2)/2) + p)*v - u*tau_xy - v*tau_yy + qy
F = np.array([F1,F2,F3,F5])
return F
#######主程序： MacCormack 方法 ############################
# 先计算下述变量，以启动循环
e = cv*T
mu = mu0*(T/T0)**1.5*(T0 + 110)/(T + 110) # 萨瑟兰公式
k = mu*cp/Pr
p = rho*R*T
U = calculate_U(rho, u, v, e)
for i in range(max_Iteration):
## 保存所需原变量
u_history[i,:,:] = u
v_history[i,:,:] = v
T_history[i,:,:] = T
rho_history[i,:,:] = rho
if i%200==0:
print('Iteration = ',i)
if (np.abs(rho-rho_history[i-1,:,:])1:#第一次迭代完密度不变
break
## 计算时间步长
nu_prime = 4.0/3.0*(gamma*mu/Pr)/rho
dt_CFL = (np.abs(u)/dx+np.abs(v)/dy+np.sqrt(gamma*R*T)*np.sqrt(1.0/dx**2 + 1.0/dy**2)+2*nu_prime*(1.0/dx**2 + 1.0/dy**2))**(-1)
dt = np.min(K*dt_CFL[1:-1,1:-1]) # 只考虑内部网格点
###### 预估步 ##################################
tau_xx_1 = tau_xx(mu,u,v,case=1) # 1 代表向后差分
tau_yy_1 = tau_yy(mu,u,v,case=1)
tau_xy_E_1 = tau_xy_E(mu,u,v,case=1)
tau_xy_F_1 = tau_xy_F(mu,u,v,case=1)
qx_1 = qx(T,k,case=1)
qy_1 = qy(T,k,case=1)
E = calculate_E(rho, u, v, e, p, tau_xx_1, tau_xy_E_1, qx_1)
F = calculate_F(rho, u, v, e, p, tau_yy_1, tau_xy_F_1, qy_1)
U_pred[:,:-1,:-1] = U[:,:-1,:-1] - dt/dx*(E[:,1:,:-1]-E[:,:-1,:-1]) - dt/dy*(F[:,:-1,1:]-F[:,:-1,:-1])
# 原变量的预测值
rho_pred,u_pred,v_pred,T_pred = calculate_original(U_pred)
# 利用边界条件完善预测值
u_pred[:,-1] = u_inf # 上边界
v_pred[:,-1] = 0
T_pred[:,-1] = T_inf
rho_pred[:,-1] = rho_inf
u_pred[-1,1:-1] = 2*u_pred[-2,1:-1]-u_pred[-3,1:-1] # 出流边界
v_pred[-1,1:-1] = 2*v_pred[-2,1:-1]-v_pred[-3,1:-1]
T_pred[-1,1:-1] = 2*T_pred[-2,1:-1]-T_pred[-3,1:-1]
rho_pred[-1,1:-1] = 2*rho_pred[-2,1:-1]-rho_pred[-3,1:-1]
u_pred[-1,0] = 0
v_pred[-1,0] = 0
T_pred[-1,0] = Tw
# T_pred[-1,0] = (4*T_pred[-1,1]-T_pred[-1,2])/3 # 绝热壁
rho_pred[-1,0] = 2*rho_pred[-1,1]-rho_pred[-1,-2]
# 计算其他变量
e_pred = cv*T_pred
mu_pred = mu0*(T_pred/T0)**1.5*(T0 + 110)/(T_pred + 110) # 萨瑟兰公式
k_pred = mu_pred*cp/Pr
p_pred = rho_pred*R*T_pred
###### 校正步 ##############################################
tau_xx_0 = tau_xx(mu_pred,u_pred,v_pred,case=0) # 0 代表向前差分
tau_yy_0 = tau_yy(mu_pred,u_pred,v_pred,case=0)
tau_xy_E_0 = tau_xy_E(mu_pred,u_pred,v_pred,case=0)
tau_xy_F_0 = tau_xy_F(mu_pred,u_pred,v_pred,case=0)
qx_0 = qx(T_pred,k_pred,case=0)
qy_0 = qy(T_pred,k_pred,case=0)
E_pred = calculate_E(rho_pred, u_pred, v_pred, e_pred, p_pred, tau_xx_0, tau_xy_E_0, qx_0)
F_pred = calculate_F(rho_pred, u_pred, v_pred, e_pred, p_pred, tau_yy_0, tau_xy_F_0, qy_0)
# 这里只算内部网格点
U[:,1:-1,1:-1] = 0.5*(U[:,1:-1,1:-1] + U_pred[:,1:-1,1:-1]
- dt/dx*(E_pred[:,1:-1,1:-1]-E_pred[:,:-2,1:-1])
- dt/dy*(F_pred[:,1:-1,1:-1]-F_pred[:,1:-1,:-2]))
rho,u,v,T = calculate_original(U)
# 边界条件，入流边界和上边界保持不变
u[-1,1:-1] = 2*u[-2,1:-1]-u[-3,1:-1] # 出流边界
v[-1,1:-1] = 2*v[-2,1:-1]-v[-3,1:-1]
T[-1,1:-1] = 2*T[-2,1:-1]-T[-3,1:-1]
rho[-1,1:-1] = 2*rho[-2,1:-1]-rho[-3,1:-1]
#下边界
T[1:,0] = Tw
# T[1:,0] = (4*T[1:,1]-T[1:,2])/3 # 绝热壁,二阶精度
rho[1:,0] = 2*rho[1:,1]-rho[1:,-2]
# 计算其他变量
e = cv*T
mu = mu0*(T/T0)**1.5*(T0 + 110)/(T + 110) # 萨瑟兰公式
k = mu*cp/Pr
p = rho*R*T
# 计算边界 U, 入流边界和上边界保持不变
U[0,:,0] = rho[:,0]
U[1,:,0] = rho[:,0]*u[:,0]
U[2,:,0] = rho[:,0]*v[:,0]
U[3,:,0] = rho[:,0]*(e[:,0]+(u[:,0]**2+v[:,0]**2)/2)
U[0,-1,1:-1] = rho[-1,1:-1]
U[1,-1,1:-1] = rho[-1,1:-1]*u[-1,1:-1]
U[2,-1,1:-1] = rho[-1,1:-1]*v[-1,1:-1]
U[3,-1,1:-1] = rho[-1,1:-1]*(e[-1,1:-1]+(u[-1,1:-1]**2+v[-1,1:-1]**2)/2)
# 去掉多余的第一维度
u_history = u_history[:i+1,:,:]
v_history = v_history[:i+1,:,:]
T_history = T_history[:i+1,:,:]
rho_history = rho_history[:i+1,:,:]
## 验证质量守恒
m_in = dy*np.sum(rho[0,1:]*u[0,1:])
m_out = dy*np.sum(rho[-1,1:]*u[-1,1:])
if np.abs(m_in-m_out) |