【游戏流体力学基础及Unity代码（十三）】压力泊松方程,python写SIMPLE算法

快乐人L · 发表于 2021-1-31 17:46

SIMPLE算法全称为Semi-Implicit Method For Pessure-Linked Equation，即求解压力耦合方程组的半隐式算法。要在python正确实现该算法，就需要知道雅可比迭代法以及压力泊松方程。本章所有代码都在下面仓库里
交错网格

交错网格(Staggered Grid)就是指下面这种东西，流体所在的格子为蓝色，假如X方向n个格子，Y方向m个格子，那么总个数为nm。流体所在格子外再围一圈用于计算压力，有(n+2)(m+2)。但速度不存在格子上，而是在压力格子之间，X轴方向的速度有(n+1)(m+2)个，Y轴方向有(n+2)(m+1)。
把速度放在格子之间这种做法就是很符合直觉，算起来也很方便。对于一维交错网格，压力网格编号与速度网格编号如下
一维泊松压力方程

泊松压力方程(Poisson Pressure Equation)在网格法中很重要，因为它能帮助准确算出没有散度的速度场，符合不可压缩流体的连续性方程。除此之外，泊松方程在上一章的涡方法和之后的ParticleInCell方法中也会用到。很多教材上会写的数学公式如下，其中尖头向下的三角形是算梯度/散度的，尖头向上的是拉普拉斯算子

$\nabla^2p = \Delta p = \nabla \cdot U\\$
一维情况下的偏微分方程和离散化如下

$\frac{\partial^2 p}{\partial x^2} = \frac{\partial u}{\partial x} \qquad \frac{p_{i+1} - 2p_{i} + p_{i-1}}{(\Delta x)^2} = \frac{u_{i+1/2} - u_{i-1/2}}{\Delta x}\\$
最标准的一维泊松压力方程应该是这样的
import numpy as np
nmax = 7
tmax = 1000
v = np.zeros((nmax-1))
v2 = np.zeros((nmax<span class="o">-1))
div = np.zeros((nmax))
div2 = np.zeros((nmax))
dx = 0.1
for i in range(1,nmax-2):
v = i+1#这步假设了一个初始有散度的速度场
for i in range(1,nmax-1):
div = (v-v[i-1])/dx
p = np.zeros((tmax,nmax))
for t in range(0,tmax-2):
for i in range(1,nmax-1):
      term = (p[t,i-1] + p[t,i+1] - dx*dx*div)/2
      p[t+1,i] = p[t,i] + (term - p[t,i])
for i in range(0,nmax-1):
v2 = v + (p[t,i] - p[t,i+1])/dx
for i in range(1,nmax-1):#重新计算散度，应当是零
div2 = (v2-v2[i-1])/dx所谓压力泊松方程就是，第一步计算一个有散度的速度场的散度，第二步将散度代入泊松方程计算，重复第二步许多次，计算出最后的压力场。第三步将初始速度场减去这个压力场的梯度，就得到了没有散度的速度场，这个速度场就符合不可压缩流体的连续性方程。
如果能正确算出压力的话，最后各点的速度v2会是一样的，这样才能保证一维情况下的没有散度的条件。最后div2算出来也是0，精度误差会让零变成大概是10的负8~17次方。
然后就是方程本身，写成代码有三种方法，一种是雅可比(jacobi)迭代，是要求计算出上一步所有变量后，才计算下一步所有变量。另一种是高斯赛德尔(Gauss-Seidel)迭代，计算出一个变量后，立马把它当成解并计算下一个变量。后者比前者速度稍微快一点。
#第一种方法，使用临时变量
for i in range(1,nmax-1)
p_star[t+1,i] = (p[t,i-1] + p[t,i+1] - div)/2
p = p_star.copy()
#第二种方法，不使用临时变量
p[t+1,i] = (p[t,i-1] + p[t,i+1] - div)/2第三种方法是求矩阵，矩阵大概看起来如下，这种矩阵又被称为三对角矩阵(Tridiagonal matrix)。不用循环，只需求一次逆矩阵即可。但为了快速求出逆矩阵还需要各种其它技巧，如The tridiagonal matrix algorithm (TDMA)。

$\frac{1}{(\Delta x)^2}\begin{bmatrix} -2 & 1 & 0 & 1 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 1 & 0 & 1 & -2\end{bmatrix}\begin{bmatrix}p_1 \\ p_2 \\ p_3 \\ p_4 \end{bmatrix} = \begin{bmatrix}div_1 \\ div_2 \\ div_3 \\ div_4 \end{bmatrix}\\$
如果不用压力的梯度为零的边界条件，而让边界压力一直为固定值的话，那么迭代很多次后压力能够收敛。但是如果应用了梯度为零的边界条件，我发现要压力最后收敛的一个条件是，在压力迭代之前的总散度要是零。例如下面的代码，1+2-3=0。压力能够正确收敛。但如果更改div[4] = -2，让总散度不是零，压力就无法收敛越算越大。我没有在任何资料上看到过这里的解释，所以我也不知道这东西叫什么。希望有知道的大佬告诉我。
import numpy as np
nmax = 7
tmax = 1000
v = np.zeros((nmax-1))
div = np.zeros((nmax))
p = np.zeros((tmax,nmax))
div[2] = 1
div[3] = 2
div[4] = -3
for t in range(0,tmax-2):
for i in range(1,nmax-1):
      term = (p[t,i-1] + p[t,i+1] - div)/2
      p[t+1,i] = p[t,i] + (term - p[t,i])
p[t+1,0] = p[t+1,1]
p[t+1,nmax-1] = p[t+1,nmax-2]这还没完，将泊松压力方程稍微变化一下，即可求解一般的线性方程组。大意就是每个网格之间的压力差权重不一样。主要代码如下，
left = np.ones((nmax))
right = np.ones((nmax))
left[4] = 2
right[3] = 2
for t in range(0,tmax-2):
for i in range(1,nmax-1):
      term = (left*p[t,i-1] + right*p[t,i+1] - dx*dx*div)/(left + right)
      p[t+1,i] = p[t,i] + (term - p[t,i])
for i in range(0,nmax-1):
v2 = v - right*(p[t,i+1] - p[t,i])*dx
for i in range(1,nmax-1):#重新计算散度，应当接近零
div2 = (v2-v2[i-1])/dx
#div2 = div - (p[t,i-1] - p[t,i])*left/dx - (p[t,i+1] - p[t,i])*right/dx完整代码在https://github.com/clatterrr/FluidSimulationTutorialsUnity/blob/main/python/SIMPLE/PoissonWeight1d.py。上面代码解的一般的线性方程组如下

$(p_0 - p_1) + (p_2 - p_1) = 1\\ (p_1 - p_2) + (p_3 - p_2) = 1\\ (p_2 - p_3) + 2(p_4 - p_3) = 1\\ 2(p_3 - p_4) + (p_5 - p_4) = 1\\ (p_4 - p_5) + (p_6 - p_5) = -4\\ 固定条件 p_0 = p_6 = 0 \\$ 为了验证正确性，我们来手算一遍上面的代码。程序计算出下面最终的压力值
散度值和之前一样没改变
由于压力高的地方的水流会流向压力低的地方，那么对于第一个网格，会向第零个网格流(-1.545 - 0) = -1.545，向第二个网格流(-1.545 -( -2.09)) = 0.545，这样一来，第一个网格的散度增加了(-1.545 + 0.545) = -1.00刚好把第一个网格的散度抵消了。由于第一个网格流向第零个网格-1.545个压力，导致这两个压力网格中间的速度网格也增加1.545个速度，注意第零个网格的压力要比第一个网格的压力高，所以速度肯定是正的。假设这里密度是一。最后的无散速度场如下
而对于第三个和第四个压力网格，它们之间的流量权重经过修改已经变成了2，也就是left[4] = right[3] = 2。如果不这样对称，压力最后虽然也能收敛，但是速度无论怎样都算不出来。在对称的情况下，第三个网格压力网格会流向第四个压力网格(-1.63 - (-0.90))*2 = -1.5个压力，而这两个网格之间的第三个速度网格的速度原本是3，加上-1.5个压力后就是1.5个速度了。其它地方也一样。最后得到无散速度场，每处的速度都是1.54，每处的散度都是零。
还有个细节是，如果一开始散度绝对值越大，那么最后算出的压力绝对值也越大。而如果压力差权重越大，那么最后算出的压力绝对值越小。你可以试着把在迭代算压力前把散度乘上十倍，看看最后压力会发生什么变化。这点在调试SIMPLE算法时也很重要。
由于压力差才是速度改变的原因，因此大多数时候压力的绝对值意义并不大，只需关注压力差就行了。
二维泊松压力方程

二维交错网格编号如下
二维泊松压力方程的偏微分形式如下

$\frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2}= \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\\$
写成离散化如下

$\frac{p_{i+1,j} - 2p_{i,j} + p_{i-1,j}}{(\Delta x)^2} + \frac{p_{i,j+1} - 2p_{i,j} + p_{i,j-1}}{(\Delta y)^2}= div\\$
再更换顺序

$p_{i,j} = \frac{(p_{i+1,j} + p_{i-1,j})(\Delta y)^2 + (p_{i,j+1} + p_{i,j-1})(\Delta x)^2- div(\Delta x)^2(\Delta y)^2}{2((\Delta x)^2 + (\Delta y)^2)}\\$
与一维相比就是加了个维度而已，其余没什么特别的。主要部分如下
for i in range(1,nmax+1):
for j in range(1,mmax+1):
      div[i,j] = (u[i,j] - u[i-1,j])/dx + (v[i,j] - v[i,j-1])/dy
for t in range(0,tmax-2):
for i in range(1,nmax+1):
      for j in range(1,mmax+1):
         term = (p[i+1,j] + p[i-1,j])*dy*dy + (p[i,j+1] + p[i,j-1])*dx*dx
         term2 = dx*dx*dy*dy*div[i,j]
         p[i,j] = p[i,j] + ((term - term2)/(2*(dx*dx + dy*dy)) - p[i,j])
for i in range(0,nmax+1):
for j in range(0,mmax+2):
      u2[i,j] = u[i,j] + (p[i,j] - p[i+1,j])/dx
for i in range(0,nmax+2):
for j in range(0,mmax+1):
      v2[i,j] = v[i,j] + (p[i,j] - p[i,j+1])/dy运行poisson2d.py代码，同样可以手算结果验证正确性。这里红字1是下面，3是上面，2是左边，4是右边。dx = 0.25,dy = 0.1666，那么对于压力网格[1,1]，这里的散度原本是10。现在往下散度增加(-0.25 - 0)/(dydy) = -9.32，往上(-0.25 - (-0.34))/(dydy) = 3.23，往左就是-0.25/(dxdx) = -4.14，往右约是0.23，如此一来正好将散度减为零。
接下来给这些压力差加上权重，这也是SIMPLE算法的核心部分。首先写离散化的泊松方程，下式的e,w,n,s分别是东方，西方，北方，南方的意思。

$\frac{e(p_{i+1,j} - p_{i,j}) + w(p_{i-1,j}- p_{i,j})}{(\Delta x)^2} + \frac{n(p_{i,j+1} - p_{i,j}) + s(p_{i,j-1}- p_{i,j})}{(\Delta y)^2}= div\\$
经过变换最终如下

$p_{i,j} = \frac{(ep_{i+1,j} + wp_{i-1,j})(\Delta y)^2 + (np_{i,j+1} + sp_{i,j-1})(\Delta x)^2- div(\Delta x)^2(\Delta y)^2}{(n + s)(\Delta x)^2 + (e + w)(\Delta y)^2}$
与之前的无权重的相比改动如下，完整代码在https://github.com/clatterrr/FluidSimulationTutorialsUnity/blob/main/python/SIMPLE/PoissonWeight2d.py
for t in range(0,tmax-2):
for i in range(1,nmax+1):
      for j in range(1,mmax+1):
         term = (e[i,j]*p[i+1,j] + w[i,j]*p[i-1,j])*dy*dy
         term2 = (n[i,j]*p[i,j+1] + s[i,j]*p[i,j-1])*dx*dx
         term3 = (e[i,j] + w[i,j])*dy*dy + (n[i,j] + s[i,j])*dx*dx
         p[i,j] = (term + term2 - dx*dx*dy*dy*div[i,j])/term3
for i in range(0,nmax+1):
for j in range(0,mmax+2):
      u2[i,j] = u[i,j] + (p[i,j] - p[i+1,j])/dx*e[i,j]
for i in range(0,nmax+2):
for j in range(0,mmax+1):
      v2[i,j] = v[i,j] + (p[i,j] - p[i,j+1])/dy*n[i,j]最后总结一下，要迭代计算压力，一般用jacobi迭代，或gauss-seidel迭代，要保证能够算出正确压力的条件是压力差权重是对称的。如果不追求物理真实，可以不加边界条件。
边界条件

在考虑物理真实的情况下，在墙壁处，压力的边界条件一般设为梯度为零，这样才能保证墙壁法向量上流体速度为零从而不穿墙。例如对于下图左，深蓝色为墙上的压力，未算出，而浅蓝色为内部压力，已算出。应用梯度为零的条件，则深蓝色网格的值也是5，很简单。
对于二维顶盖驱动流的边界，比如区域的左下角如上图中，三个深蓝色网格都填5，也能满足梯度为零条件。但对于方形绕流，比如方形障碍物的右上角如上图右，深蓝色网格的压力值应该填什么？
我搜了一圈资料都没找到这东西应该怎么填。接下来的方法是我自己想出来的，可能与别人有出入。但是目前来看效果不错。我的方法就是，如果要处理如果左侧是边界压力，就暂时让左侧压力梯度为零。如果要处理下侧边界压力，就暂时让下侧压力梯度为零。主要代码如下，其中block=1的区域为障碍物。
for i in range(1,nmax+1):
      for j in range(1,mmax+1):
         if(block[i,j] == 1):
            continue
         if(block[i-1,j] == 1):
            p[i-1,j] = p[i,j]
         if(block[i,j-1] == 1):
            p[i,j-1] = p[i,j]
         term = (p[i+1,j] + p[i-1,j])*dy*dy + (p[i,j+1] + p[i,j-1])*dx*dx
         term -= dx*dx*dy*dy*div[i,j]
         p[i,j] = p[i,j] + (term/(2*(dx*dx + dy*dy)) - p[i,j])完整代码在
SIMPLE算法

关于SIMPLE算法的很多参考资料都只写了公式没代码，尽管如此，关于SIMPLE算法的介绍不算少，因此接下来仅仅展示如何把SIMPLE算法写成代码。SIMPLE算法可分为如下四步
第一步，猜测速度场u,v，猜一个压力场p
第二步，用u,v解速度动量方程，得到u_star,v_star，这个速度场有散度
第三步，用u_star,v_star,p解泊松压力方程，得到p_star
第四步，用p_star计算得到正确的速度和压力场，这个速度场无散度
第二步首先是离散化动量方程，同样SIMPLE算法用ewns表示东西北南方向，我也不知道为什么有这样的习惯。

$\frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u^2)}{\partial x} = -\frac{\partial p}{\partial x} + \frac{1}{Re}(\frac{\partial^2 u}{\partial x^2})\\ \frac{\partial (\rho u)}{\partial t} = -((\rho u)_eu_e - (\rho u)_wu_w) + \frac{1}{Re}((\frac{\partial u}{\partial x})_e - (\frac{\partial u}{\partial x})_w)-\frac{\partial p}{\partial x}\\$
然后将参数提取出来。至于为什么这么写各位还是看参考文献去吧，这玩意没几千字解释不清楚...

$A_e = \frac{1}{Re} - \frac{F_e}{2} \quad A_w = \frac{1}{Re} + \frac{F_e}{2}\quad A_n = \frac{1}{Re} - \frac{F_n}{2} \quad A_s = \frac{1}{Re} + \frac{F_s}{2}\\ Fe = max(0,-\rho u_e) = \frac{|\rho u_e| - \rho u_e}{2} \\ Fw = max(0,\rho u_w) = \frac{|\rho u_w| + \rho u_w}{2}\\$
这样动量方程就变成了如下形式，或者称为Rhie-Chow插值

$A_pu_p = A_eu_e + A_wu_w + \nabla p = \sum_{ew}Au + \nabla p\\$
其中

$A_p = \frac{1}{Re} + \frac{1}{Re} + \frac{F_e}{2} - \frac{F_w}{2} = A_e + A_w + F_e - F_w\\$
二维动量方程经过上面的方法同样可以变成下面这样子，其实与之前说的泊松压力方程非常类似。计算中心点的速度用到了东南西北四个方向，所以被称为半隐式的。

$A_pu_p = A_eu_e + A_wu_w + A_nu_n + A_su_s + \nabla p = \sum_{ewns}Au + \Delta y(p_e - p_2)\\$
有些参考资料上还会写成这种形式

$u_p = \frac{\sum_{ewns}Au}{A_p} + d(p_e - p_2) \qquad d = \frac{\Delta y}{A_p}\\$
这些AeAwAnAs也可以看成是速度差权重，比如对于u11，位置如下
上面形式的动量方程写成代码如下
for i in range(1, nmax):
for j in range(1, mmax + 1):
      term = Ae[i, j] * u_star[i + 1, j] + Aw[i, j] * u_star[i - 1, j
      <span class="n">term += As[i, j] * u_star[i, j - 1] + An[i, j] * u_star[i, j + 1]
      pressure_term = (p[i, j] - p[i + 1, j]) * dy
      u_star[i, j] = u_star[i, j] + a*((term + pressure_term) / Apu[i, j] - u_star[i, j])其中参数计算方法
Fe = 0.5 * (u[i, j] + u[i + 1, j]) * dy * rho
Fw = 0.5 * (u[i, j] + u[i - 1, j]) * dy * rho
Fn = 0.5 * (v[i, j] + v[i + 1, j]) * dx * rho
Fs = 0.5 * (v[i, j - 1] + v[i + 1, j - 1]) * dx * rho
Ae[i, j] = max(0, -Fe) + dy / dx / Re
Aw[i, j] = max(0, Fw) + dy / dx / Re
An[i, j] = max(0, -Fn) + dx / dy / Re
As[i, j] = max(0, Fs) + dx / dy / Re
Apu[i, j] = Ae[i, j] + Aw[i, j] + An[i, j] + As[i, j] + Fe - Fw + Fn - Fs计算v_star方法是一样的。然后是正牌的压力校正法，和之前的离散动量方程，以及泊松压力方程差不多。在很多SIMPLE资料中它长这样

$ap = a_ep_e + a_wp_w + a_np_n + a_sp_s - b_{i,j}= \sum_{ewns} ap - b_{i,j}\\$
写成代码如下，相比于动量方程就是改了个变量名字
term = Ape[i, j] * p_star[i + 1, j] + Apw[i, j] * p_star[i - 1, j]
term += Apn[i, j] * p_star[i, j + 1] + Aps[i, j] * p_star[i, j - 1]
p_source[i, j] =(u_star[i, j]-u_star[i - 1, j])*dy + (v_star[i, j] - v_star[i, j - 1]) * dx
term = (term - p_source[i, j]) / App[i, j] # p_source就是公式那个b，也就是速度的散度
p_star[i, j] = p_star[i, j] + 0.8 * (term - p_star[i, j])其中的系数计算如下，或者说是压力差的权重计算如下。其具体含义也请读者查阅相关资料。

$a_e = \rho d_e \Delta y \quad a_w =\rho d_w \Delta y \quad a_n =\rho d_n \Delta x \quad a_s =\rho d_s \Delta x\\ d_e = \frac{\Delta y}{A_{pe}} \quad d_w = \frac{\Delta y}{A_{pw}} \quad d_n = \frac{\Delta x}{A_{pn}} \quad d_w = \frac{\Delta x}{A_{ps}}\\$
写成代码如下
for i in range(1, nmax + 1):
  for j in range(1, mmax + 1):
   if i < nmax:
      Ape[i, j] = rho * dy * dy / Apu[i, j]
   if i > 1:
      Apw[i, j] = rho * dy * dy / Apu[i - 1, j]
   if j < mmax:
      Apn[i, j] = rho * dx * dx / Apv[i, j]
   if j > 1:
      Aps[i, j] = rho * dx * dx / Apv[i, j - 1]
   App[i, j] = Ape[i, j] + Apw[i, j] + Apn[i, j] + Aps[i, j]这样算法主要部分就完了。不过还有个任务，请读者调试时观察Ae,Aw,An,As,Ap,Ape,Apw,Apn,Aps,App这几个数组里数据的特点，然后解释它为什么会出现这样的特点。SIMPLE解决顶盖驱动流问题的完整代码在下
https://github.com/clatterrr/FluidSimulationTutorialsUnity/blob/main/python/SIMPLE/SIMPLEcavity.py
如果是管道流，与顶盖驱动流不同是初始速度条件。
# Pipe Flow初始速度条件，保证入口总速度和出口总速度一样
u[0, 1:mmax+1] = 1
u[nmax, 2:mmax] = mmax/(mmax-2)  以及每个时间步中压力的边界条件，进出口的压力是固定值，具体取什么值要看实际的物理情况。最后速度结果如下。因为入口固定为6个1，出口固定为中间4个1.5。
后台阶流动也只需要修改一下初始速度条件和压力边界条件。最后部分速度如下，下图红圈处就是后台阶的位置，后台阶的前面有个顺时针的漩涡。
然后是方块绕流，处理压力边界条件之前已经说过。本来是想模拟卡门涡街的，结果并没模拟出，原因暂时没找到。最终结果如下，左侧X轴方形速度，右侧y轴方形速度。中间那一堆速度为零的地方就是方形障碍物的位置。
总结

SIMPLE算法我调试代码一个月才慢慢理解了。虽然SIMPLE主要是CFD里用，但自己写一遍也能学到不少东西。我调试的主要方法如下
1.首先去github上把与此有关的代码都下载下来，然后去搜索相关的文章，尽快把公式和代码给联系上。不仅可以在github网站内部搜索，也可以直接在外部搜索引擎搜索“xxx github”,搜索“xxx matlab”有可能在mathwork网站上找到代码。另外一些代码可能在cfdonline/油管视频/某些附代码的书中。尽量用不同的关键字试着搜索，比如有些SIMPLE算法的资料是我找PISO算法的资料时才找出来的。
2.如果看文章看不懂的话，有人选择把它抄一遍。我觉得这并不是好方法。我一般会找二三十篇资料都看一遍，然后选出五六篇易懂的继续反复看，再看不懂就继续找新资料。因为有的资料中忽略的地方在别的资料中可能说的很详细。
3.自己写的时候，先从简单的写起，比如先一维后二维再三维，先4x4网格或4个粒子，然后再增加网格粒子数量。然后验证简单的算例，比如顶盖流，管道流，后台阶流，圆柱绕流。
4.去理解每个参数的物理意义，思考每一步的结果是否既物理正确也数值正确，比如压力泊松和SIMPLE算法每个时间步结束后，都要保证速度的散度处处为零
5.程序算一遍，然后自己手算一遍验算，防止程序漏掉某项。我写泊松压力方程时一开始漏掉个除dx，手算了一下午才找出来
6.如果程序写得正确，那么它是从第一个时间步开始就是正确的。比如用SIMPLE算法模拟顶盖驱动流，只需要一个时间步，就可以算出合理的速度场和压力场。
7.收敛。我学了SIMPLE算法才发现泊松压力方程的重要性，为了研究这个方程什么条件下压力才会收敛而不会越算越大，我把这个方程单独提取了出来重新写。于是就有了本篇的第二节和第三节。
8.多写代码。SIMPLE我只找到了几份顶盖流的代码可参考，资料算是不少。以后某些算法资料稀少，参考代码约等于没有，这时候就得全凭自己过去的思考了。
9.还是写不出来干脆不写了。等过几个月或几年经验积累足够了，或是突然找到一份绝佳的资料，会发现以前的难题其实非常简单
本篇文章没写到的地方包括
1.SIMPLE算法不仅可以在交错网格上，也可以放在规则网格上(Collected Grid)上，后者速度定义在网格中间。
2.如果你觉得SIMPLE算法太简单，可以挑战一下困难模式，实现SIMPLER,SIMPLEC,PISO,PIMPLE算法，参考资料更少，而且至少我没找到参考代码
接下来是可压缩流体的内容，可压缩流体在超声速飞机与火箭喷管中都很重要。
上一篇：光影帽子：【游戏流体力学基础及Unity代码（十二）】卡门涡街，边界层，涡方法
参考

《热流体数值计算方法与应用》作者：何志霞，王谦，袁建平，第二章解线性方程组，第五章SIMPLE算法写得非常详细。附有c语言代码
http://dyfluid.com/simplefoam.html OpenFoam中的simpleFoam模块解析
《Numerical Heat Transfer and Fluid Flow》 by Suhas V Patankar 请认准这本书的作者，他就是SIMPLE算法的发明者。此书第六章介绍了SIMPLE和SIMPLER，写得很详细
《The Finite Volume Method in Computational Fluid Dynamics An Advanced Introduction with OpenFOAM and Matlab》 by F. Moukalled, L. Mangani, M. Darwish第15章，主要说了交错网格和非交错网格上的SIMPLE，以及SIMPLEC和PISO，有一些手算例子
《Computational Fluid Mechanics and Heat Transfer》by Anderson, Dale Pletcher, Richard H. Tannehill, John C 第九章，介绍了SIMPLE,PISO以及在交错网格上的Collected Grid上的写法
《数值方法(matlab版，第四版)》,第三章的雅可比迭代介绍得很好
Pressure Correction Scheme for Incompressible Fluid Flow  by Ong Chin Kai
Solution of Navier-Stokes Equations for Incompressible Flows Using SIMPLE and MAC Algorithms
https://github.com/davidlin001/flow_sims
https://www.mathworks.com/matlabcentral/fileexchange/68348-2d-lid-driven-cavity-flow-using-simple-algorithm
https:/github.com/deepmorzaria/Lid-Driven-Cavity-Collocated-Grid
https://github.com/Vignesh2508/NS_Solver_SIMPLE
https://github.com/Gnakman/Python_Numerics

		自动登录	找回密码
密码			立即注册

[笔记] 【游戏流体力学基础及Unity代码（十三）】压力泊松方程,python写SIMPLE算法

本帖子中包含更多资源