【游戏流体力学基础及Unity代码（十三）】压力泊松方程,python写SIMPLE算法
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方法中也会用到。很多教材上会写的数学公式如下，其中尖头向下的三角形是算梯度/散度的，尖头向上的是拉普拉斯算子
https://www.zhihu.com/equation?tex=%5Cnabla%5E2p+%3D+%5CDelta+p+%3D+%5Cnabla+%5Ccdot+U%5C%5C
一维情况下的偏微分方程和离散化如下
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial%5E2+p%7D%7B%5Cpartial+x%5E2%7D+%3D+%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D+%5Cqquad+%5Cfrac%7Bp_%7Bi%2B1%7D+-+2p_%7Bi%7D+%2B+p_%7Bi-1%7D%7D%7B%28%5CDelta+x%29%5E2%7D+%3D+%5Cfrac%7Bu_%7Bi%2B1%2F2%7D+-+u_%7Bi-1%2F2%7D%7D%7B%5CDelta+x%7D%5C%5C
最标准的一维泊松压力方程应该是这样的
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)/dx
p = np.zeros((tmax,nmax))
for t in range(0,tmax-2):
for i in range(1,nmax-1):
term = (p + p - dx*dx*div)/2
p = p + (term - p)
for i in range(0,nmax-1):
v2 = v + (p - p)/dx
for i in range(1,nmax-1):#重新计算散度，应当是零
div2 = (v2-v2)/dx所谓压力泊松方程就是，第一步计算一个有散度的速度场的散度，第二步将散度代入泊松方程计算，重复第二步许多次，计算出最后的压力场。第三步将初始速度场减去这个压力场的梯度，就得到了没有散度的速度场，这个速度场就符合不可压缩流体的连续性方程。
如果能正确算出压力的话，最后各点的速度v2会是一样的，这样才能保证一维情况下的没有散度的条件。最后div2算出来也是0，精度误差会让零变成大概是10的负8~17次方。
然后就是方程本身，写成代码有三种方法，一种是雅可比(jacobi)迭代，是要求计算出上一步所有变量后，才计算下一步所有变量。另一种是高斯赛德尔(Gauss-Seidel)迭代，计算出一个变量后，立马把它当成解并计算下一个变量。后者比前者速度稍微快一点。
#第一种方法，使用临时变量
for i in range(1,nmax-1)
p_star = (p + p - div)/2
p = p_star.copy()
#第二种方法，不使用临时变量
p = (p + p - div)/2第三种方法是求矩阵，矩阵大概看起来如下，这种矩阵又被称为三对角矩阵(Tridiagonal matrix)。不用循环，只需求一次逆矩阵即可。但为了快速求出逆矩阵还需要各种其它技巧，如The tridiagonal matrix algorithm (TDMA)。
https://www.zhihu.com/equation?tex=%5Cfrac%7B1%7D%7B%28%5CDelta+x%29%5E2%7D%5Cbegin%7Bbmatrix%7D+-2+%26+1+%26+0+%26+1+%5C%5C+1+%26+-2+%26+1+%26+0+%5C%5C+0+%26+1+%26+-2+%26+1+%5C%5C+1+%26+0+%26+1+%26+-2%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7Dp_1+%5C%5C+p_2+%5C%5C+p_3+%5C%5C+p_4++%5Cend%7Bbmatrix%7D+%3D+%5Cbegin%7Bbmatrix%7Ddiv_1+%5C%5C+div_2+%5C%5C+div_3+%5C%5C+div_4++%5Cend%7Bbmatrix%7D%5C%5C
如果不用压力的梯度为零的边界条件，而让边界压力一直为固定值的话，那么迭代很多次后压力能够收敛。但是如果应用了梯度为零的边界条件，我发现要压力最后收敛的一个条件是，在压力迭代之前的总散度要是零。例如下面的代码，1+2-3=0。压力能够正确收敛。但如果更改div = -2，让总散度不是零，压力就无法收敛越算越大。我没有在任何资料上看到过这里的解释，所以我也不知道这东西叫什么。希望有知道的大佬告诉我。
import numpy as np
nmax = 7
tmax = 1000
v = np.zeros((nmax-1))
div = np.zeros((nmax))
p = np.zeros((tmax,nmax))
div = 1
div = 2
div = -3
for t in range(0,tmax-2):
for i in range(1,nmax-1):
term = (p + p - div)/2
p = p + (term - p)
p = p
p = p这还没完，将泊松压力方程稍微变化一下，即可求解一般的线性方程组。大意就是每个网格之间的压力差权重不一样。主要代码如下，
left = np.ones((nmax))
right = np.ones((nmax))
left = 2
right = 2
for t in range(0,tmax-2):
for i in range(1,nmax-1):
term = (left*p + right*p - dx*dx*div)/(left + right)
p = p + (term - p)
for i in range(0,nmax-1):
v2 = v - right*(p - p)*dx
for i in range(1,nmax-1):#重新计算散度，应当接近零
div2 = (v2-v2)/dx
#div2 = div - (p - p)*left/dx - (p - p)*right/dx完整代码在https://github.com/clatterrr/FluidSimulationTutorialsUnity/blob/main/python/SIMPLE/PoissonWeight1d.py。上面代码解的一般的线性方程组如下
https://www.zhihu.com/equation?tex=%28p_0+-+p_1%29+%2B+%28p_2+-+p_1%29+%3D+1%5C%5C+%28p_1+-+p_2%29+%2B+%28p_3+-+p_2%29+%3D+1%5C%5C+%28p_2+-+p_3%29+%2B+2%28p_4+-+p_3%29+%3D+1%5C%5C+2%28p_3+-+p_4%29+%2B+%28p_5+-+p_4%29+%3D+1%5C%5C+%28p_4+-+p_5%29+%2B+%28p_6+-+p_5%29+%3D+-4%5C%5C+%E5%9B%BA%E5%AE%9A%E6%9D%A1%E4%BB%B6+p_0+%3D+p_6+%3D+0+%5C%5C 为了验证正确性，我们来手算一遍上面的代码。程序计算出下面最终的压力值
散度值和之前一样没改变
由于压力高的地方的水流会流向压力低的地方，那么对于第一个网格，会向第零个网格流(-1.545 - 0) = -1.545，向第二个网格流(-1.545 -( -2.09)) = 0.545，这样一来，第一个网格的散度增加了(-1.545 + 0.545) = -1.00刚好把第一个网格的散度抵消了。由于第一个网格流向第零个网格-1.545个压力，导致这两个压力网格中间的速度网格也增加1.545个速度，注意第零个网格的压力要比第一个网格的压力高，所以速度肯定是正的。假设这里密度是一。最后的无散速度场如下
而对于第三个和第四个压力网格，它们之间的流量权重经过修改已经变成了2，也就是left = right = 2。如果不这样对称，压力最后虽然也能收敛，但是速度无论怎样都算不出来。在对称的情况下，第三个网格压力网格会流向第四个压力网格(-1.63 - (-0.90))*2 = -1.5个压力，而这两个网格之间的第三个速度网格的速度原本是3，加上-1.5个压力后就是1.5个速度了。其它地方也一样。最后得到无散速度场，每处的速度都是1.54，每处的散度都是零。
还有个细节是，如果一开始散度绝对值越大，那么最后算出的压力绝对值也越大。而如果压力差权重越大，那么最后算出的压力绝对值越小。你可以试着把在迭代算压力前把散度乘上十倍，看看最后压力会发生什么变化。这点在调试SIMPLE算法时也很重要。
由于压力差才是速度改变的原因，因此大多数时候压力的绝对值意义并不大，只需关注压力差就行了。
二维泊松压力方程
二维交错网格编号如下
二维泊松压力方程的偏微分形式如下
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial%5E2+p%7D%7B%5Cpartial+x%5E2%7D+%2B+%5Cfrac%7B%5Cpartial%5E2+p%7D%7B%5Cpartial+y%5E2%7D%3D+%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D+%2B+%5Cfrac%7B%5Cpartial+v%7D%7B%5Cpartial+y%7D%5C%5C
写成离散化如下
https://www.zhihu.com/equation?tex=%5Cfrac%7Bp_%7Bi%2B1%2Cj%7D+-+2p_%7Bi%2Cj%7D+%2B+p_%7Bi-1%2Cj%7D%7D%7B%28%5CDelta+x%29%5E2%7D+%2B+%5Cfrac%7Bp_%7Bi%2Cj%2B1%7D+-+2p_%7Bi%2Cj%7D+%2B+p_%7Bi%2Cj-1%7D%7D%7B%28%5CDelta+y%29%5E2%7D%3D+div%5C%5C
再更换顺序
https://www.zhihu.com/equation?tex=p_%7Bi%2Cj%7D+%3D+%5Cfrac%7B%28p_%7Bi%2B1%2Cj%7D++%2B+p_%7Bi-1%2Cj%7D%29%28%5CDelta+y%29%5E2+%2B+%28p_%7Bi%2Cj%2B1%7D+%2B+p_%7Bi%2Cj-1%7D%29%28%5CDelta+x%29%5E2-+div%28%5CDelta+x%29%5E2%28%5CDelta+y%29%5E2%7D%7B2%28%28%5CDelta+x%29%5E2+%2B+%28%5CDelta+y%29%5E2%29%7D%5C%5C
与一维相比就是加了个维度而已，其余没什么特别的。主要部分如下
for i in range(1,nmax+1):
for j in range(1,mmax+1):
div = (u - u)/dx + (v - v)/dy
for t in range(0,tmax-2):
for i in range(1,nmax+1):
for j in range(1,mmax+1):
term = (p + p)*dy*dy + (p + p)*dx*dx
term2 = dx*dx*dy*dy*div
p = p + ((term - term2)/(2*(dx*dx + dy*dy)) - p)
for i in range(0,nmax+1):
for j in range(0,mmax+2):
u2 = u + (p - p)/dx
for i in range(0,nmax+2):
for j in range(0,mmax+1):
v2 = v + (p - p)/dy运行poisson2d.py代码，同样可以手算结果验证正确性。这里红字1是下面，3是上面，2是左边，4是右边。dx = 0.25,dy = 0.1666，那么对于压力网格，这里的散度原本是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分别是东方，西方，北方，南方的意思。
https://www.zhihu.com/equation?tex=%5Cfrac%7Be%28p_%7Bi%2B1%2Cj%7D+-+p_%7Bi%2Cj%7D%29+%2B+w%28p_%7Bi-1%2Cj%7D-+p_%7Bi%2Cj%7D%29%7D%7B%28%5CDelta+x%29%5E2%7D+%2B+%5Cfrac%7Bn%28p_%7Bi%2Cj%2B1%7D+-+p_%7Bi%2Cj%7D%29+%2B+s%28p_%7Bi%2Cj-1%7D-+p_%7Bi%2Cj%7D%29%7D%7B%28%5CDelta+y%29%5E2%7D%3D+div%5C%5C
经过变换最终如下
https://www.zhihu.com/equation?tex=p_%7Bi%2Cj%7D+%3D+%5Cfrac%7B%28ep_%7Bi%2B1%2Cj%7D++%2B+wp_%7Bi-1%2Cj%7D%29%28%5CDelta+y%29%5E2+%2B+%28np_%7Bi%2Cj%2B1%7D+%2B+sp_%7Bi%2Cj-1%7D%29%28%5CDelta+x%29%5E2-+div%28%5CDelta+x%29%5E2%28%5CDelta+y%29%5E2%7D%7B%28n+%2B+s%29%28%5CDelta+x%29%5E2+%2B+%28e+%2B+w%29%28%5CDelta+y%29%5E2%7D
与之前的无权重的相比改动如下，完整代码在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*p + w*p)*dy*dy
term2 = (n*p + s*p)*dx*dx
term3 = (e + w)*dy*dy + (n + s)*dx*dx
p = (term + term2 - dx*dx*dy*dy*div)/term3
for i in range(0,nmax+1):
for j in range(0,mmax+2):
u2 = u + (p - p)/dx*e
for i in range(0,nmax+2):
for j in range(0,mmax+1):
v2 = v + (p - p)/dy*n最后总结一下，要迭代计算压力，一般用jacobi迭代，或gauss-seidel迭代，要保证能够算出正确压力的条件是压力差权重是对称的。如果不追求物理真实，可以不加边界条件。
边界条件
在考虑物理真实的情况下，在墙壁处，压力的边界条件一般设为梯度为零，这样才能保证墙壁法向量上流体速度为零从而不穿墙。例如对于下图左，深蓝色为墙上的压力，未算出，而浅蓝色为内部压力，已算出。应用梯度为零的条件，则深蓝色网格的值也是5，很简单。
对于二维顶盖驱动流的边界，比如区域的左下角如上图中，三个深蓝色网格都填5，也能满足梯度为零条件。但对于方形绕流，比如方形障碍物的右上角如上图右，深蓝色网格的压力值应该填什么？
我搜了一圈资料都没找到这东西应该怎么填。接下来的方法是我自己想出来的，可能与别人有出入。但是目前来看效果不错。我的方法就是，如果要处理如果左侧是边界压力，就暂时让左侧压力梯度为零。如果要处理下侧边界压力，就暂时让下侧压力梯度为零。主要代码如下，其中block=1的区域为障碍物。
for i in range(1,nmax+1):
for j in range(1,mmax+1):
if(block == 1):
continue
if(block == 1):
p = p
if(block == 1):
p = p
term = (p + p)*dy*dy + (p + p)*dx*dx
term -= dx*dx*dy*dy*div
p = p + (term/(2*(dx*dx + dy*dy)) - p)完整代码在
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表示东西北南方向，我也不知道为什么有这样的习惯。
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+%28%5Crho+u%29%7D%7B%5Cpartial+t%7D+%2B+%5Cfrac%7B%5Cpartial+%28%5Crho+u%5E2%29%7D%7B%5Cpartial+x%7D+%3D+-%5Cfrac%7B%5Cpartial+p%7D%7B%5Cpartial+x%7D+%2B+%5Cfrac%7B1%7D%7BRe%7D%28%5Cfrac%7B%5Cpartial%5E2+u%7D%7B%5Cpartial+x%5E2%7D%29%5C%5C+%5Cfrac%7B%5Cpartial+%28%5Crho+u%29%7D%7B%5Cpartial+t%7D++%3D+-%28%28%5Crho+u%29_eu_e+-+%28%5Crho+u%29_wu_w%29+%2B+%5Cfrac%7B1%7D%7BRe%7D%28%28%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D%29_e+-+%28%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D%29_w%29-%5Cfrac%7B%5Cpartial+p%7D%7B%5Cpartial+x%7D%5C%5C
然后将参数提取出来。至于为什么这么写各位还是看参考文献去吧，这玩意没几千字解释不清楚...
https://www.zhihu.com/equation?tex=A_e+%3D+%5Cfrac%7B1%7D%7BRe%7D+-+%5Cfrac%7BF_e%7D%7B2%7D+%5Cquad+A_w+%3D+%5Cfrac%7B1%7D%7BRe%7D+%2B+%5Cfrac%7BF_e%7D%7B2%7D%5Cquad++A_n+%3D+%5Cfrac%7B1%7D%7BRe%7D+-+%5Cfrac%7BF_n%7D%7B2%7D+%5Cquad+A_s+%3D+%5Cfrac%7B1%7D%7BRe%7D+%2B+%5Cfrac%7BF_s%7D%7B2%7D%5C%5C+Fe+%3D+max%280%2C-%5Crho+u_e%29+%3D+%5Cfrac%7B%7C%5Crho+u_e%7C+-+%5Crho+u_e%7D%7B2%7D+%5C%5C++Fw+%3D+max%280%2C%5Crho+u_w%29+%3D+%5Cfrac%7B%7C%5Crho+u_w%7C+%2B+%5Crho+u_w%7D%7B2%7D%5C%5C
这样动量方程就变成了如下形式，或者称为Rhie-Chow插值
https://www.zhihu.com/equation?tex=A_pu_p+%3D+A_eu_e+%2B+A_wu_w+%2B+%5Cnabla+p+%3D+%5Csum_%7Bew%7DAu+%2B+%5Cnabla+p%5C%5C
其中
https://www.zhihu.com/equation?tex=A_p+%3D+%5Cfrac%7B1%7D%7BRe%7D+%2B+%5Cfrac%7B1%7D%7BRe%7D+%2B+%5Cfrac%7BF_e%7D%7B2%7D+-+%5Cfrac%7BF_w%7D%7B2%7D+%3D+A_e+%2B+A_w+%2B++F_e+-+F_w%5C%5C
二维动量方程经过上面的方法同样可以变成下面这样子，其实与之前说的泊松压力方程非常类似。计算中心点的速度用到了东南西北四个方向，所以被称为半隐式的。
https://www.zhihu.com/equation?tex=A_pu_p+%3D+A_eu_e+%2B+A_wu_w+%2B+A_nu_n+%2B+A_su_s+%2B+%5Cnabla+p+%3D+%5Csum_%7Bewns%7DAu+%2B+%5CDelta+y%28p_e+-+p_2%29%5C%5C
有些参考资料上还会写成这种形式
https://www.zhihu.com/equation?tex=u_p+%3D++%5Cfrac%7B%5Csum_%7Bewns%7DAu%7D%7BA_p%7D+%2B+d%28p_e+-+p_2%29+%5Cqquad+d+%3D+%5Cfrac%7B%5CDelta+y%7D%7BA_p%7D%5C%5C
这些AeAwAnAs也可以看成是速度差权重，比如对于u11，位置如下
上面形式的动量方程写成代码如下
for i in range(1, nmax):
for j in range(1, mmax + 1):
term = Ae * u_star + Aw * u_star[i - 1, j
<span class="n">term += As * u_star + An * u_star
pressure_term = (p - p) * dy
u_star = u_star + a*((term + pressure_term) / Apu - u_star)其中参数计算方法
Fe = 0.5 * (u + u) * dy * rho
Fw = 0.5 * (u + u) * dy * rho
Fn = 0.5 * (v + v) * dx * rho
Fs = 0.5 * (v + v) * dx * rho
Ae = max(0, -Fe) + dy / dx / Re
Aw = max(0, Fw) + dy / dx / Re
An = max(0, -Fn) + dx / dy / Re
As = max(0, Fs) + dx / dy / Re
Apu = Ae + Aw + An + As + Fe - Fw + Fn - Fs计算v_star方法是一样的。然后是正牌的压力校正法，和之前的离散动量方程，以及泊松压力方程差不多。在很多SIMPLE资料中它长这样
https://www.zhihu.com/equation?tex=ap+%3D+a_ep_e+%2B+a_wp_w+%2B+a_np_n+%2B+a_sp_s+-+b_%7Bi%2Cj%7D%3D+%5Csum_%7Bewns%7D+ap+-+b_%7Bi%2Cj%7D%5C%5C
写成代码如下，相比于动量方程就是改了个变量名字
term = Ape * p_star + Apw * p_star
term += Apn * p_star + Aps * p_star
p_source =(u_star-u_star)*dy + (v_star - v_star) * dx
term = (term - p_source) / App # p_source就是公式那个b，也就是速度的散度
p_star = p_star + 0.8 * (term - p_star)其中的系数计算如下，或者说是压力差的权重计算如下。其具体含义也请读者查阅相关资料。
https://www.zhihu.com/equation?tex=a_e+%3D+%5Crho+d_e+%5CDelta+y+%5Cquad+a_w+%3D%5Crho+d_w+%5CDelta+y+%5Cquad+a_n+%3D%5Crho+d_n+%5CDelta+x+%5Cquad+a_s+%3D%5Crho+d_s+%5CDelta+x%5C%5C+d_e+%3D+%5Cfrac%7B%5CDelta+y%7D%7BA_%7Bpe%7D%7D+%5Cquad+d_w+%3D+%5Cfrac%7B%5CDelta+y%7D%7BA_%7Bpw%7D%7D+%5Cquad+d_n+%3D+%5Cfrac%7B%5CDelta+x%7D%7BA_%7Bpn%7D%7D+%5Cquad+d_w+%3D+%5Cfrac%7B%5CDelta+x%7D%7BA_%7Bps%7D%7D%5C%5C
写成代码如下
for i in range(1, nmax + 1):
for j in range(1, mmax + 1):
if i < nmax:
Ape = rho * dy * dy / Apu
if i > 1:
Apw = rho * dy * dy / Apu
if j < mmax:
Apn = rho * dx * dx / Apv
if j > 1:
Aps = rho * dx * dx / Apv
App = Ape + Apw + Apn + Aps这样算法主要部分就完了。不过还有个任务，请读者调试时观察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 = 1
u = 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 Flowby 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
页:
[1]