不可压库艾特流的压力修正法求解（附完整代码）

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

实现了 第 9.4 节另一种数值方法：压力修正法 的代码。该方法通过采用三套不同的网格，分别在不同的网格上计算压力，速度水平分量，速度垂直分量，从而避免出现所求物理量出现“棋盘式分布”。因为采用了交错网格，记录起来比较麻烦，所以采用了Python中的3阶ndarray来表示网格组，具体说明如下。

1、网格的示意图如下所示，其中黑点代表“压力(p)网格”，红叉代表“垂直速度(v)网格”，蓝点代表“水平速度(u)网格”。以黑点为基准对整个计算域进行网格划分，而后红叉在垂直方向上与黑点交错，蓝点在水平方向上与黑点交错。如图所示，用一个直角三角形将（黑点，蓝点，红叉）三个构成一组，最外一圈用铅笔实线勾勒出的网格点组表示计算边界，通过边界条件可得，因此计算时，只需根据公式计算内部的网格点组，在Python中，x方向和y方向对应的数组切片均是[1:-1]。

2、具体而言，构造了一个3阶的 ndarray ，大小为 (xnum,ynum,3)，第一个参数表示 x 方向的网格组数，第二参数表示 y 方向的网格组数，第三个参数表示要求解的参数，分别是 0: p，1: u ，2: v。如果要记录迭代过程中的值，可以再增加一个维度，表示迭代次数。采用此种方式表示，相应的公式下标较原书中有所变化（书中采用 1/2 表示，其实并不方便编程），现将差分方程中要用到的式子重新表述如下：

$\frac{\partial (\rho u^2)}{\partial x} = \frac{\rho u_{i+1,j}^2 - \rho u_{i-1,j}^2}{2\Delta x}$ ； $\frac{\partial (\rho v^2)}{\partial y} = \frac{\rho v_{i,j+1}^2 - \rho v_{i,j-1}^2}{2\Delta y}$

$\frac{\partial (\rho u v)}{\partial y} = \frac{\rho u_{i+1,j} \bar{v} - \rho u_{i-1,j} v}{2\Delta y}$ ，其中 $\bar{v} = \frac{v_{i,j} + v_{i+1,j}}{2}$ ， $v = \frac{v_{i,j-1} + v_{i+1,j-1}}{2}$

$\frac{\partial (\rho u v)}{\partial x} = \frac{\rho v_{i+1,j} \bar{u} - \rho v_{i-1,j} u}{2\Delta x}$ ，其中 $\bar{u} = \frac{u_{i,j} + u_{i,j+1}}{2}$ ， $u = \frac{u_{i-1,j} +u_{i-1,j+1}}{2}$

$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{(\Delta x)^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{(\Delta y)^2}$

$\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} = \frac{v_{i+1,j} - 2v_{i,j} + v_{i-1,j}}{(\Delta x)^2} + \frac{v_{i,j+1} - 2v_{i,j} + v_{i,j-1}}{(\Delta y)^2}$

$\frac{\partial p}{\partial x} = \frac{p_{i+1,j} - p_{i,j} }{\Delta x}$ ； $\frac{\partial p}{\partial y} = \frac{p_{i,j+1} - p_{i,j} }{\Delta y}$

压力修正公式中质量源项 d 的计算式：

$d = \frac{\rho u_{i,j}^* - \rho u_{i-1,j}^* }{\Delta x} + \frac{\rho v_{i,j} - \rho v_{i,j-1}^* }{\Delta y}$

代码运行结果如下：

完整代码如下：

# -*- coding: utf-8 -*-
"""
Created on Fri Aug 14 00:41:09 2020@author: PANSONG
"""##############################################
######## 不可压库艾特流动：压力修正法 ##########
##############################################import numpy as np
import matplotlib.pyplot as plt;plt.close('all')#%% 几何条件
L = 0.5 # ft，x 方向长
D = 0.01 # ft, y 方向高### 物性
rho = 0.002377 # slug/ft3 密度
mu = 3.74e-7 ## ReD = 63.6 粘度### 边界条件
ue = 1 # ft/s 上边界速度
ve = 0 # 下边界速度
pe_prime = 0 # 上边界压力修正量u0 = 0 # 下边界条件
v0 = 0
p0_prime = 0p_inlet_prime = 0 # 入流边界
v_inlet = 0p_outlet_prime = 0 # 出流边界### 生成网格
xnum = 21 # x 方向网格数，以 p 为基准
ynum = 11 # y 方向网格数dx = L/(xnum-1) # x 方向网格宽度
dy = D/(ynum-1) # y 方向网格宽度
dt = 0.001y = np.linspace(0,D,ynum)pset = np.zeros((xnum,ynum,3)) # 用三阶数组表示，前两个参数表示x,y方向，3代表3个参数 0:p,1:u,2:v;### 计算参数,及临时变量
max_iteration = 300 # 最大主迭代次数
max_sub_iteration = 500 # 最大子迭代次数p_prime = np.zeros((xnum,ynum)) # p'
p_prime_old = p_prime.copy() # 用于判断松弛法计算 p' 是否收敛pset_history = np.zeros((max_iteration,xnum,ynum,3)) # 保存 p,u,v 迭代值
d_history = np.zeros(max_iteration)#%% 计算：压力修正法，SIMPLE 算法 ###################
## 第一步：初始化
pset[:,-1,1] = ue
pset[:,-1,2] = ve
pset[:,0,1] = u0
pset[:,0,2] = v0
pset[0,:,2] = v_inletpset[15,5,2] = ue/2 # 点（15,5）的 v 脉冲p_prime[0,:] = p_inlet_prime
p_prime[-1,:] = p_outlet_prime
p_prime[:,0] = p0_prime
p_prime[:,-1] = pe_primefor i in range(max_iteration):pset_history[i,:,:,:] = pset # 保存第 i 次迭代值## 第二步：计算 rhou*,rhov*v_bar = 0.5*(pset[1:-1,1:-1,2] + pset[2:,1:-1,2])v = 0.5*(pset[1:-1,0:-2,2] + pset[2:,0:-2,2])A_star = ( -(rho*pset[2:,1:-1,1]**2 - rho*pset[:-2,1:-1,1]**2)/dx/2.0 -(rho*pset[1:-1,2:,1]*v_bar - rho*pset[1:-1,:-2,1]*v)/dy/2.0 + mu*(pset[2:,1:-1,1] - 2*pset[1:-1,1:-1,1] + pset[:-2,1:-1,1])/(dx**2)+ mu*(pset[1:-1,2:,1] - 2*pset[1:-1,1:-1,1] + pset[1:-1,:-2,1])/(dy**2) )u_bar= 0.5*(pset[1:-1,1:-1,1] + pset[1:-1,2:,1])u = 0.5*(pset[0:-2,1:-1,1] + pset[0:-2,2:,1])B_star = ( -(rho*pset[1:-1,2:,2]**2 - rho*pset[1:-1,:-2,2]**2)/dy/2.0 -(rho*pset[2:,1:-1,2]*u_bar - rho*pset[:-2,1:-1,2]*u)/dx/2.0 + mu*(pset[2:,1:-1,2] - 2*pset[1:-1,1:-1,2] + pset[:-2,1:-1,2])/(dx**2)+ mu*(pset[1:-1,2:,2] - 2*pset[1:-1,1:-1,2] + pset[1:-1,:-2,2])/(dy**2) )pset[1:-1,1:-1,1] = 1.0/rho*(rho*pset[1:-1,1:-1,1] + A_star*dt - dt/dx*(pset[2:,1:-1,0] - pset[1:-1,1:-1,0]))pset[1:-1,1:-1,2] = 1.0/rho*(rho*pset[1:-1,1:-1,2] + B_star*dt - dt/dy*(pset[1:-1,2:,0] - pset[1:-1,1:-1,0]))## 零阶外插，得到入流边界和出流边界值;上下边界保持不变， v_inlet恒为0pset[0,1:-1,1] = pset[1,1:-1,1] pset[-1,1:-1,1] = pset[-2,1:-1,1]pset[-1,1:-1,2] = pset[-2,1:-1,2]    ## 第三步：计算压力修正量a = 2*(dt/(dx**2) + dt/(dy**2))b = - dt/(dx**2)c = - dt/(dy**2)d = (rho*pset[1:-1,1:-1,1]-rho*pset[:-2,1:-1,1])/dx + (rho*pset[1:-1,1:-1,2]-rho*pset[1:-1,:-2,2])/dyd_history[i] = d[14,4] # 保存点[15,5]处质量源项的变化for j in range(max_sub_iteration):if j == (max_sub_iteration-1):print('Warning ! max iteration for p_prime')p_prime[1:-1,1:-1] = -1.0/a*(b*p_prime[2:,1:-1]+b*p_prime[:-2,1:-1]+c*p_prime[1:-1,2:]+c*p_prime[1:-1,:-2]+d) # 更新 p', p' 的边界值保持不变if np.linalg.norm((p_prime-p_prime_old))<1e-9: #用二阶范数判断,相当于距离breakp_prime_old = p_prime.copy()## 第四步：更新 p 值pset[:,:,0] = pset[:,:,0] + 0.1*p_prime # 低松弛因子取 0.1#%% 结果可视化 #######
## 网格点 i=15 处，速度分量 v 随垂直高度 y的分布
K_list = [0,1,4,50,100,-1] # 不同迭代次数
K_legend = []
for K in K_list:plt.plot(pset_history[K,15,:,2],y)K_legend.append('K='+str(K))plt.legend(K_legend)
plt.xlabel('v (ft/s)')
plt.ylabel('y (ft)')
plt.title('Evolution of v in the middle')## 质量源项 d
plt.figure()
plt.plot(d_history)
plt.xlabel('Iteration')
plt.ylabel('d (slug/ft$^3$·s)')
plt.title('Evolution of mass source in (15,5)')## i = 15, u 的变化
plt.figure()
K_list2 = [4,20,50,150,299]
K_legend2 = []
for K in K_list2:plt.plot(pset_history[K,15,:,1],y)K_legend2.append('K='+str(K))plt.legend(K_legend2)
plt.xlabel('u (ft/s)')
plt.ylabel('y (ft)')
plt.title('Evoluion of u in the middle')