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

入门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 = 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')

 

 

  • 5
    点赞
  • 19
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
以下是使用Python实现不可压库艾特流完整代码: ```python import numpy as np import matplotlib.pyplot as plt # 定义流体的参数 rho = 1.225 # 空气密度 U = 10 # 流体速度 # 定义物体的参数 c = 1 # 翼弦长 alpha = 5 # 迎角 N = 100 # 离散化的分段数 # 将物体离散化为N个小分段 dz = c/N # 小分段的长度 z = dz/2 + np.arange(0, N)*dz # 定义边界条件 u_inf = U*np.ones(N) v_inf = np.zeros(N) # 定义迭代参数 tol = 1e-6 # 迭代收敛的精度 maxiter = 100 # 最大迭代次数 # 初始化迭代变量 u = u_inf v = v_inf err = 1 # 迭代求解 for i in range(maxiter): u_old = u.copy() v_old = v.copy() # 计算涡度 omega = np.zeros(N) for j in range(N): omega[j] = (u[(j+1)%N] - u[j])/dz - (v[(j+1)%N] + v[j])/c # 计算压力系数 cp = 1 - (u**2 + v**2)/U**2 # 计算总的升力系数 cl = np.dot(cp, np.sin(alpha*np.pi/180)*dz/c) # 根据涡量方程求解速度场 for j in range(N): u[j] = u_inf[j] + 1/(2*np.pi)*sum((v[(j-k-1)%N] + v[(j-k)%N])/2*dz/np.sqrt((z[j]-z[k])**2+c**2) for k in range(N)) v[j] = v_inf[j] - 1/(2*np.pi)*sum((u[(j-k-1)%N] + u[(j-k)%N])/2*dz/np.sqrt((z[j]-z[k])**2+c**2) for k in range(N)) + omega[j]*c/2 # 判断迭代是否收敛 err = np.max(np.abs(u - u_old)) if err < tol: break # 输出结果 print('Lift coefficient: ', cl) # 绘制压力系数分布图 plt.plot(z/c, cp) plt.xlabel('z/c') plt.ylabel('Cp') plt.gca().invert_yaxis() plt.show() ``` 在这个示例中,我们使用迭代方求解不可压库塔-祖考夫斯基定理。首先,我们将物体离散化为N个小分段,并定义流体和物体的参数。然后,我们根据涡量方程求解速度场,并计算压力系数和总的升力系数。最后,我们绘制压力系数分布图并输出结果。
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值