入门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 表示,其实并不方便编程),现将差分方程中要用到的式子重新表述如下:
;
,其中
,
,其中
,
;
压力修正公式中质量源项 d 的计算式:
代码运行结果如下:
完整代码如下:
# -*- 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 = 0
p_inlet_prime = 0 # 入流边界
v_inlet = 0
p_outlet_prime = 0 # 出流边界
### 生成网格
xnum = 21 # x 方向网格数,以 p 为基准
ynum = 11 # y 方向网格数
dx = L/(xnum-1) # x 方向网格宽度
dy = D/(ynum-1) # y 方向网格宽度
dt = 0.001
y = 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_inlet
pset[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_prime
for 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恒为0
pset[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])/dy
d_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: #用二阶范数判断,相当于距离
break
p_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')