SIMPLE算法求解多孔介质的一维流动控制方程

问题介绍

设流经某多孔介质的一维流动控制方程为：
$$C|u|u+dp/dx=0$$$$d(uF)/dx=0$$
其中，系数$C$与空间位置有关，$F$为流道的有效截面积。
对于图1所示的均匀网格，假定：
$$C_B=0.25$$$$C_C=0.2$$$$F_B=5$$$$F_C=4$$$$p_1=200$$$$p_3=38$$$$\Delta x=2$$
以上各单位都是协调的（也就是说不用管单位，直接带入数值计算就得了），试采用$SIMPLE$算法确定$p_2、u_B、u_C$的值？

之前对于$|u|u$的写法很不理解，20210803突然想到，这样的写法不就体现出方向性了吗？如果写成$u^2$就不能体现方向性了。所以这样的写法本人猜测其核心是为了体现方向性。

求解思想

首先，从图1我们可以看出，本算例采用的是交错网格（staggered grid）：速度$u$和压力$p$分别存储于两套不同的网格系统；速度$u$存储于压力$p$控制容积的东、西界面上。
其次，我们主要是体会$SIMPLE$算法的精髓：压力修正+预估校正。
最后，在程序中体现迭代校正。自拟初始估计值。

压力修正方法的基本思想

1、假设一个压力场，记为$p^{\ast }$;
2、利用$p^{\ast }$，求解动量离散方程，得到相应的速度$u^{\ast }$;
3、利用质量守恒方程来改进压力场，要求与改进后的压力场相对应的速度场能满足连续方程；
4、用$p^{{}^{\prime }}$表示压力修正值，$u^{{}^{\prime }}$表示速度修正值，以$p^{\ast }+p^{{}^{\prime }}$及$u^{\ast }+u^{{}^{\prime }}$作为本层次的解，并据此开始下一层次的迭代计算。

两个关键问题

1、如何获得压力修正值$p^{{}^{\prime }}$，使（$p^{\ast }+p^{{}^{\prime }}$）满足连续性方程?
2、获得$p^{{}^{\prime }}$后，如何确定$u^{{}^{\prime }}$?

求解步骤及说明

求解步骤与压力修正方法的基本思想不同，一个是先假定速度，一个是先假定压力，哪个是对的？
答：没有矛盾。先假定速度是为了求动量离散方程的系数及源项

疑惑：$a_e$表示动量方程离散系数？怎么求解？

本案例求解分析过程

本案例计算Python代码

# 已知参数
CB = 0.25
CC = 0.2
FB = 5
FC = 4
p1 = 200
p3 = 38
dx = 2

# B、C处的速度试算值和初始压力假定
uB_0 = 8
uC_0 = 10
p2_star = 100

# 设置统计迭代次数变量
count = 1

# 迭代开始
while 1:
    # 求速度预估值
    uB_star = 0.5 * uB_0 - (p2_star - p1) / (2 * CB * uB_0 * dx)
    uC_star = 0.5 * uC_0 - (p3 - p2_star) / (2 * CC * uC_0 * dx)

    # 求压力修正值
    p2_correction = 2 * dx * (uB_star * FB - uC_star * FC) / ((FB / (CB * uB_0)) + (FC / (CC * uC_0)))

    # 求速度修正值
    uB_correction = - p2_correction / (2 * CB * uB_0 * dx)
    uC_correction = p2_correction / (2 * CC * uC_0 * dx)

    # 迭代结果
    p2 = p2_star + p2_correction
    uB = uB_star + uB_correction
    uC = uC_star + uC_correction
    
    # 打印第count次的迭代结果（第count次的真值，可能不是最终真正的真值）
    print('第',count,'次迭代结果：')
    print('p2 = ',p2) 
    print('uB = ',uB)
    print('uC = ',uC)
    print('\n')
    # 退出循环判断
   # if abs(p2_correction) < 0.000001 and abs(uB_correction) < 0.000001 and abs(uC_correction) <0.000001:#退出循环条件
    if abs(p2 - p2_star) < 0.000001 and abs(uB - uB_0) < 0.000001 and abs(uC - uC_0) <0.000001:
        # 条件判断中的p2_star、uB_0、uC_0是上一迭代层的真值（从下面else中赋值语句可以看出）
        # 条件判断中的p2、uB、uC是最新迭代层（本迭代层）的真值
        # 以上一迭代层的真值与本迭代层的真值的差值小于10的负6次方作为收敛条件（收敛后每一迭代层的真值理应相等）
        break
    else:
        p2_star = p2
        uB_0 = uB
        uC_0 = uC
        count += 1
        
print('只需迭代',count,'次使得残差值都小于10的负6次方')
print('最终求得的结果为：')
print('p2 = ',p2) 
print('uB = ',uB)
print('uC = ',uC)
  

计算结果如下：