入门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 = .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+(u2+v2)/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,:,:]-(u2+v2)/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*u2 + p - tau_xx E3 = rho*u*v - tau_xy E5 = (rho*(e+(u2+v2)/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*v2 + p - tau_yy F5 = (rho*(e+(u2+v2)/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,:,:])<1e-8).all() and i>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/dx2 + 1.0/dy2)+2*nu_prime*(1.0/dx2 + 1.0/dy2))(-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)<0.01: print('') print(' 质量守恒定律满足 ！') print('') else: print('#') print(' 质量守恒定律不满足 ！') print('#') #%% 计算结果可视化 分别注释 绝热壁条件 和 等温壁条件 后运行代码，并将稳态结果保存 等温壁 # rho_Tw = rho # T_Tw = T # u_Tw = u # v_Tw = v 绝热壁 # rho_ad = rho # T_ad = T # u_ad = u # v_ad = v p_Tw = rho_Tw*R*T_Tw p_ad = rho_ad*R*T_ad # # 沿平板表面的网格站位，网格点 j=1 plt.figure() plt.plot(p_Tw[:,1]/p_inf) plt.plot(p_ad[:,1]/p_inf) plt.xlabel('Mesh position') plt.ylabel('p/$p_\infty$') plt.legend(['Isothermal','Adiabatic']) plt.title('Pressure distribution on the slab surface') # # 平板后缘压力剖面 # # 边界层坐标 Rex = rho_inf*u_inf*x[-1]/mu0 y_bar = y/x[-1]*np.sqrt(Rex) # 压力 plt.figure() plt.plot(p_Tw[-1,:]/p_inf,y_bar) plt.plot(p_ad[-1,:]/p_inf,y_bar) plt.xlabel('p/$p_\infty$ on the rear') plt.ylabel('boundary layer y coordinate') plt.legend(['Isothermal','Adiabatic']) plt.title('Pressure distribution on the slab rear') # 密度 plt.figure() plt.plot(rho_Tw[-1,:]/rho_inf,y_bar) plt.plot(rho_ad[-1,:]/rho_inf,y_bar) plt.xlabel('$\\rho$/$\\rho_\\infty$ on the rear') plt.ylabel('boundary layer y coordinate') plt.legend(['Isothermal','Adiabatic']) plt.title('Density distribution on the slab rear') # 温度 plt.figure() plt.plot(T_Tw[-1,:]/T_inf,y_bar) plt.plot(T_ad[-1,:]/T_inf,y_bar) plt.xlabel('T/$T_\infty$ on the rear') plt.ylabel('boundary layer y coordinate') plt.legend(['Isothermal','Adiabatic']) plt.title('Temperature distribution on the slab rear') # 水平速度 u plt.figure() plt.plot(u_Tw[-1,:]/u_inf,y_bar) plt.plot(u_ad[-1,:]/u_inf,y_bar) plt.xlabel('u/$u_\infty$ on the rear') plt.ylabel('boundary layer y coordinate') plt.legend(['Isothermal','Adiabatic']) plt.title('Horizontal velocity distribution on the slab rear') # 马赫数 Ma_Tw = np.sqrt(u_Tw2+v_Tw2)/np.sqrt(gamma*R*T_Tw) Ma_ad = np.sqrt(u_ad2+v_ad2)/np.sqrt(gamma*R*T_ad) plt.figure() plt.plot(Ma_Tw[-1,:],y_bar) plt.plot(Ma_ad[-1,:],y_bar) plt.xlabel('Local Mach number on the rear') plt.ylabel('boundary layer y coordinate') plt.legend(['Isothermal','Adiabatic']) plt.title('Mach number distribution on the slab rear') 等温壁条件等值图 X,Y = np.meshgrid(x,y) # 压力 plt.figure() plt.contourf(X,Y,p_Tw.T,cmap='jet') plt.xlabel('x') plt.ylabel('y') plt.title('Pressure distribution \n under isothermal wall') plt.colorbar() # 温度 plt.figure() plt.contourf(X,Y,T_Tw.T,cmap='jet') plt.xlabel('x') plt.ylabel('y') plt.title('Temperature distribution \n under isothermal wall') plt.colorbar() # 密度 plt.figure() plt.contourf(X,Y,rho_Tw.T,cmap='jet') plt.xlabel('x') plt.ylabel('y') plt.title('Density distribution \n under isothermal wall') plt.colorbar() # 马赫数 plt.figure() plt.contourf(X,Y,Ma_Tw.T,cmap='jet') plt.xlabel('x') plt.ylabel('y') plt.title('Mach number distribution \n under isothermal wall') plt.colorbar()