SIMPLE 算法基础

不可压缩流体的控制方程（直角坐标系下）：

$$ \frac{\partial \rho}{\partial t}+\frac{\partial }{\partial x_i}\left(\rho u_i\right) = 0\newline \frac{\partial \rho u_i}{\partial t}+\frac{\partial }{\partial x_i}\left(\rho u_i u_j\right) = \frac{\partial }{\partial x_i}\left(\tau_i,_j\right)-\frac{\partial p}{\partial x_i} +S_i \tag{1} $$

将方程离散后，将面临两大问题：

• 棋盘分布的流场（如下图所示）也可能满足离散后的方程组（由于压力中心格式导致）–> 采用交错网格可解决该问题（本文不讨论该话题） 。
• 耦合求解方式计算成本较大 –> 采用分离求解的方法。

SIMPLE 算法核心思想

SIMPLE 算法是主流商软的求解方法，1972年由 S. V. Patankar 和D. B. Splading 提出，算法全称为 Semi-Implicit Method for Pressure-Linked Equations（压力耦合方程组的半隐式方法），算法核心思想如下：

1. 对于给定的压力场下，可以分别求得满足压力场速度 $\bold{U}^\ast$。
2. 所得的速度场$\bold{U^\ast}$ 此时未必满足质量守恒方程，进而反向对压力场进行修正，得到 $P^1$。
3. 此时，修改后的压力场可能满足质量守恒，但同时又违反了动量守恒。因此又将 $P^1$ 带入速度修正方程，再次获得速度场 $\bold{U^1}$。
4. 返回步骤1，直至收敛。

SIMPLE 算法推导

本文整理了两种推导，其中第一种推导直观易懂，也可参考[1]。第二种较为抽象，可参考 [2]。

方式1

假设原有的速度场为 $u_0, v_0, w_0$ 和 $p_0$，。动量方程有：

$$ a_e u_e^* = \sum(a_n,_b u_n,_b,_0) + S+(p_P,_0-p_E,_0)A_e \tag{2} $$

上述方程为离散后的动量方程，详情可参阅《数值传热学（第二版）》，下面类似方程不再解释。

通过公式（2）求解的速度满足基于初始压力场 $p_0$ 的动量守恒方程，但是往往不满足连续性方程。并且由于常密度假设，连续性方程中的密度值无法更新，那么如何使得速度满足连续性方程呢？

假设即满足连续性方程又满足动量方程的物理场为： $$ u_1 = u_0 + u’ \approx u^\ast + u’,\newline v_1 = v_0 + v’ \approx v^\ast + v’,\newline w_1 = w_0 + w’ \approx w^\ast + w’,\newline p_1 = p_0 + p’ \nonumber $$ 基于 $u_1、u_1、w_1$ 和 $p_1$ 的动量方程可写作（假设体积源项不变）： $$ a_e u_e,_1 = \sum(a_n,_b u_n,_b,_1) + S+(p_P,_1-p_E,_1)A_e \tag{3} $$ 下式 (3) 减上式 (2) 可得：

$$ a_e u_e’ = \sum(a_n,_b u_n,_b’) +(p_P’-p_E’)A_e \tag{4} $$

公式（4）表明：速度的修正量 = 周边节点的速度波动 + 同向的压力波动。求解（4）的计算量较大（二维为4对角阵），且考虑到右边项的第二项 >> 右边项第一项，因此可简化为：

$$ a_e u_e’ =(p_P’-p_E’)A_e \rightarrow u_e’ =(p_P’-p_E’)(A_e/a_e)=d_e(p_P’-p_E’) \nonumber $$

类似地，可得到，

$$ v_n’ =(p_P’-p_N’)(A_e/a_e)=d_n(p_P’-p_E’) \nonumber $$

因而，修正后的速度可写为：

$$ u_e,_1=u_e^*+d_e(p_P’-p_E’),v_n,_1=u_n^\ast +d_n(p_P’-p_E’) \tag{5} $$

公式（5）称作速度的修正方程，该方程反应了速度对压力的变化。然而，物理场的连续性问题仍未解决，还是得回归连续性方程，离散后的连续性方程可写为： $$ \frac{\rho_p-\rho_P,_0}{\Delta t}+\left[\rho_e u_e -\rho_w v_w \right]\Delta y+\left[\rho_n v_n -\rho_s v_s\right]\Delta x = 0 \tag{6} $$ 为使等式成立，公式中的变量应使用 $u_1、u_1、w_1$ 和 $p_1$ 。基于上述假设，将公式（5）带入公式（6）中，并且假设初始场满足连续性方程，此时有： $$ a_Pp_P’=a_Ep_E’+a_Wp_w’+a_Np_N’+a_Sp_S’+b\ \tag{7} $$ 式中：

$$ a_E=\rho_ed_e\Delta y, \newline a_W=\rho_wd_w\Delta y, \newline a_N=\rho_nd_n\Delta x, \newline a_S=\rho_sd_s\Delta x,\newline a_P=a_E+a_W+a_N+a_S,\newline b=\frac{(\rho_P,_0 – \rho_P)\Delta x \Delta y}{\Delta t} + \left[(\rho_w u_w)^\ast-(\rho_eu_e)^\ast \right]\Delta y + \left[(\rho_s v_s)^\ast-(\rho_n v_n)^\ast\right]\Delta x. \nonumber $$

公式（7）为泊松方程，可以认为压力基于速度的变化量，或者基于原始速度的压力分布。

因此，结合上述分析，SIMPLE 算法流程为：

1. 假定一个速度分布（$\bold{U}^0$），通过计算动量方程获取方程中的系数和常数项，并假定压力场为 $p^{0}$（此时的速度场可能不满足动量方程）；
2. 分别求解各方向动量方程，获得更新的速度场 （$\bold{U}^\ast$）；
3. 基于初场（$\bold{U}^{0}$）求解压力修正方程，得到压力修正场（$p’$）；
4. 根据压力修正值，反推更新后的速度场的变化量（$\bold{U}’$）【可以理解为：根据假设和公式（5），近似出由于压力场发生变化而导致的速度场变化量】；
5. 利用改进后的速度场（$\bold{U}^1=\bold{U}^\ast+\bold{U}’$）求解那些通过源项、物性等与速度场耦合的 $\phi$ 变量。如果并不影响流场，则应在速度场收敛后再求。
6. 重新计算动量离散方程的系数，并用改进后的压力场（$p^1 = p^0+p’$）作为下一层次送代计算的初值。重复步骤2，直到获得收敛的解。

【Q&A】

1. 为什么在获得压力的改进值 $p_1=p_0+p’$后，不直接利用这一改进值及 $u^\ast$，$v^\ast$ 去开始下一层次的迭代，而是要先计算这一层次上的修正值 $u’$、$u’$呢？

   答： $u^\ast$，$u^\ast$ 不满足连续性方程，如果用它们去确定新的系数、开始下一层次的选代，则会影响选代 收敛速度，并且也使 $ a_P =\sum a_n,_b$ 的关系得不到保证，会使代数方程组系数矩阵对角占优的条件遭到破坏。[1]

   个人观点：压力改进值是基于 $u_0$、$v_0$ 进行改进的（即，$p_1$ 、 $u_0$、$v_0$ 是满足连续性方程的 ）；而 $u^\ast$，$v^\ast$ 是基于 $p_0$ 求出的（即，$p_0$ 、 $u^\ast$、$v^\ast$ 是满足动量方程的），强行将 $p_1$ 、 $u^\ast$、$v^\ast$ 可能造成更新的流场即不满足连续性方程又不满足动量方程 ，造成矩阵混乱。

方式2

假设，初场为 $\bold{U}^0$，$p^0$，修正差值为 $\bold{U}’$，$ p’$。假设动量方程的系数矩阵为$\bold{A}$ ， 基于初场 $p_0$，为满足动量方程，可求出速度场 $\bold{U}^{*}$：

$$ \bold{A^0U^*}=-\nabla p_0 \tag{8} $$

更新的速度场 $\bold{U}^\ast$ 可能不满足连续性方程，但满足动量方程。理论上，此时应该求解连续性方程，然而由于假设密度为常数，无法直接求解，因此，需要构造出新的方程使其满足连续性方程。

假设，某一流场 $f$，同时满足动量方程和连续性方程： $$ \bold{A^fU^f}=-\nabla p_0 \rightarrow \bold{D^fU^f}+ \bold{N^fU^f} = -\nabla p_f\rightarrow \bold{D^fU^f} – \bold{H}^f=-\nabla p_f \nonumber $$ 式中， $\bold{D}$ 为矩阵$\bold{A}$ 的对角阵，$\bold{N}$ 为矩阵$\bold{A}$ 的非对角阵，因此初始速度场也可以表示为： $$ \bold{U^f} =(\bold{D^f})^{-1}\bold{H}^f-(\bold{D^f})^{-1}\nabla p_f \tag{9} $$ 将等价的速度场带入连续性方程可得： $$ \nabla\cdot(\bold{D}^f)^{-1}(\nabla p_f) = \nabla\cdot(\bold{D^f})^{-1}\bold{H^f} \tag{10} $$ 公式 (10) 为满足连续性方程和动量方程的压力泊松方程。

然而，实际情况下，初始压力场不满足公式 (10)，因此，需要对初始场压力进行修正： $$ \nabla\cdot(\bold{D}^0)^{-1}(\nabla p) = \nabla\cdot(\bold{D^0})^{-1}\bold{H^0} \tag{11} $$ 式中， $$ \bold{H}^0 = (\bold{A^0 – N^0})\bold{U}^0 \nonumber $$ 根据公式 (11) 可求出压力场的变化量 $p’$，式 (11) 可以理解为满足初始速度场和连续性方程的压力场分布。

此时，求解（8）、（9）产生了一些冲突：公式 (8) 调整速度场，使其满足初始压力场；而公式 (11)，调整压力场使其满足初始速度场，但同时调整后，连续性方程和速度方程可能都不再满足。根据上述分析，要满足式 (10)，必须满足式 (9)，因此，可以根据式 (9) 求出速度的二次更新量： $$ \bold{U^1} =(\bold{D^0})^{-1}\bold{H}^0-(\bold{D^f})^{-1}\nabla p_1 \rightarrow \bold{U}’ =-(\bold{D^f})^{-1}\nabla p’ \tag{12} $$ 综上，方式2 推导的 SIMPLE 算法的流程为（其实与方式1一样）：

1. 根据动量方程求解速度场，这里的速度场并不满足连续方程：$\bold{A^0U^*}=-\nabla p_0$；
2. 根据泊松方程求解压力场：$\nabla\cdot(\bold{D}^0)^{-1}(\nabla p) = \nabla\cdot(\bold{D^0})^{-1}\bold{H^0}$；
3. 根据压力场去修正速度场，使其满足连续性方程：$\bold{U^1} =(\bold{D^\ast})^{-1}\bold{H}^\ast-(\bold{D^f})^{-1}\nabla p_1 $ 或 $\bold{U}’ =-(\bold{D^\ast})^{-1}\nabla p’$；
4. 修正后的速度场如果仍不满足连续性方程，就需要重新第一步，进行循环。

参考资料

[1] 陶文铨.数值传热学(第2版)[M].西安交通大学出版社,2001.

[2] https://www.zhihu.com/tardis/bd/art/370351186?source_id=1001

Views: 46