0%

SIMPLE，SIMPLEC和PISO算法

Posted on Edited on In

SIMPLE

算法概述

$$
{ c } \frac { ( \rho \mathbf { v } ) ^ { n + 1 } - ( \rho \mathbf { v } ) ^ { n } } { \Delta t } + C \left( \mathbf { v } ^ { n + 1 } \right) = L \left( \mathbf { v } ^ { n + 1 } \right) - G \left( p ^ { n + 1 } \right) \tag{1}
$$

$$
\frac { 3 ( \rho \mathbf { v } ) ^ { n + 1 } - 4 ( \rho \mathbf { v } ) ^ { n } + ( \rho \mathbf { v } ) ^ { n - 1 } } { 2 \Delta t } + C \left( \mathbf { v } ^ { n + 1 } \right) = L \left( \mathbf { v } ^ { n + 1 } \right) - G \left( p ^ { n + 1 } \right)\tag{2}
$$

第一个方程：隐式Euler方案；第二个方程：时间导数近似

在一个时间步长的迭代中，非线性项和耦合项被更新，称为outer iteration，以区别于在固定系数的线性系统上执行的内部迭代。

当迭代收敛时，可以得到$\mathbf v^{n+1}\approx \mathbf v^m$

连续求解线性化动量方程，当所有隐性项被组合在一起时，对于每个速度分量，我们可以得到如下的矩阵方程
$$
A^{m-1}u_i^m=Q_i^{m-1}-G_i(p^m)\tag{3}
$$

$G_i$是第i个分量的梯度；源项Q包含所有关于$u^{m - 1}_ i$可被显式求解的项，或这些项可能及任何体积力或其他线性化或延迟修正项，或在新的时间步依赖于$u^{n+1}_ i$以的其他变量（例如T）

上面方程可被改为：
$$
A_Pu_{i,P}^m+\sum_kA_ku_{i,k}^m=Q_P-(\frac{\delta p^m}{\delta x_i})_P\tag{4}
$$

压力项以符号差分形式表示，以强调解的方法与空间导数的离散化近似的独立性。
注意，系数$A_k$包含离散化对流和扩散项的贡献，而对角系数$A_P$包含非定常项的贡献

分离（3）中的A为对角部分$A_D$和非对角线部分$A_{OD}$，在外迭代m中，使用上一个时间步求得的压力
$$
(A_D+A_{OD})u_i^=Q-G_i(p^{m-1})\tag{5}
$$
通过求解（5）可求得速度场，但是不满足连续性方程，需进行速度，压力修正
$$
p^=p^{m-1}+p’,\quad u^{}_i=u_i^{*}+u_i^{’}\tag{6}
$$
采用上式获得的速度压力修正，忽略相邻点的修正，可得到
$$
A_Du_i^{}+A_{OD}u_i^{}=Q-G_i(p^)\tag{7}
$$
（7）-（5），可得
$$
A_Du_i’=-G_i(p’)\Rightarrow u_i’=-(A_D)^{-1}G_i(p’)\tag{8}
$$

  • 如果在将修正速度应用到非对角项，求解会非常的复杂

  • 上述简化会影响收敛速率，因此需要under-relaxation factors

当前，我们需要修正速度$u_i^{**}$满足离散的连续性方程，因此可得到
$$
D(\rho\mathbf v^*)+D(\rho\mathbf v’)=0\tag{9}
$$

D代表导数，不可压缩流体中$D(\rho v^m)=0$

将（8）代入（9），我们可以得到
$$
D(\rho(A_D)^{-1}G_i(p’))=D(\rho\mathbf v^*)\tag{10}
$$
上述方法通常称为投影法：首先构造一个不满足连续性方程的速度场，然后通过减去一些东西(通常是一个压力梯度)来修正它。

从向量数学的角度看，压力作为一个算子通过连续性约束将散度速度向量场投射到非发散向量场中

一旦压力修正被求解，速度和压力可通过（8）和（6）更新

上式即为**SIMPLE（Semi-Implicit Method for Pressure-Linked Equations）**算法

计算步骤

  1. 通过上一时间步m-1的速度和压力求解矩阵方程，解得矩阵系数$A_D,A_{OD}$和源项$Q^{m-1}$
  2. 通过式（5）计算得到预测速度$u_i^*$
  3. 将得到的预测速度代入泊松方程（10）中计算得到压力修正项$p’$
  4. 代入（8）式中得到速度修正项$u’_i$
  5. 将求得的速度压力修正项代入到式（6）中求得下一时间步m的速度和压力$u_i^m, p^m$
  6. 重复步骤1-5，开始迭代

SIMPLEC(SIMPLE_Corrected)

并没有忽略对于非对角项的速度修正，而是采用加权平均临近点值的方法近似速度修正
$$
u_{i,P}’\approx\frac{\sum_kA_ku_{i,k}’}{\sum_kA_k}\rightarrow \sum_kA_ku_{i,k}’\approx u_{i,P}’\sum_kA_k\tag{11}
$$
（4）减去（5）得
$$
A_Du_i’+A_{OD}u_i’=-G_i(p’)\tag{12}
$$
将（11）代入（12）中，可以得到
$$
u_i’=-(A_D+\widetilde{A_{OD}})^{-1}G_i(p’)\tag{13}
$$

$\widetilde A_D$代表对角线元素的和

ps：

  1. 对于有限容积理论，当仅存在对流项和扩散项时，存在关系式$(A_D+\widetilde{A_{OD}})=0$，然而，当存在非稳态项时，或当除以松弛因子时，都会对$A_D$产生影响，而$A_D$（通常为正）远大于$\widetilde A_{OD}$（通常为负）的绝对值
  2. 对于SIMPLEC算法，$(A_D+\widetilde{A_{OD}})$通常远小于SIMPLE算法中$A_D$的值（因为在SIMPLE算法中，压力梯度项要除以$A_D$的值，所以要避免大数除以小数。这也是SIMPLE算法中在压力修正时需要松弛因子的原因。

计算压力修正项
$$
D(\rho(A_D+\widetilde{A_{OD}})^{-1}G_i(p’))=D(\rho\mathbf v^*)\tag{14}
$$

PISO

介绍：将SIMPLE算法的步骤设定为预测步骤，同时加入一系列的修正步骤

在SIMPLE后的第一次修正，要求修正后的速度和压力满足连续性方程
$$
A_Du_i^{}+A_{OD}u_i^{}=Q-G_i(p^{})\tag{15}
$$
（15）减去（7），得到
$$
A_Du_i’’+A_{OD}u_i’=-G_i(p’)\Rightarrow u_i’’=-(A_D)^{-1}(A_{OD}u_i’+G_i(p’’))\tag{16}
$$
$u_i^{}$满足连续性方程
$$
D(\rho\mathbf v^{})+D(\rho\mathbf v’’)=0\tag{17}
$$
将（16）代入（17），同时由于$D(\rho\mathbf v^{})=0$，可得到第二个压力-修正方程**
$$
D(\rho(A_D)^{-1}G_i(p’’))=D(\rho(A_D)^{-1}A_{OD}u_i’)\tag{18}
$$
通过以上步骤重复，多次进行压力修正，即为PISO算法（Pressure Implicit with Splitting of Operators）