最近感觉CSDN有点大病，用各种浏览器都不显示最上面那一栏（就是点赞、数据、发布啥啥的那一栏），今天用ubuntu偶然发现又显示了， 赶紧把之前写的东西发出来记录一下，不知道问题出在哪：(

CFDpython - 12 steps to N-S equation

后面四步，每一步都是一个不同的方程，分别是：二维拉普拉斯方程、二维泊松方程、二维空腔流动、二维管渠流动

Step 9: 2D Laplace Equation

1. 方程形式如下：
   ∂ 2 p ∂ x 2 + ∂ 2 p ∂ y 2 = 0 \\frac{\\partial ^2 p}{\\partial x^2} + \\frac{\\partial ^2 p}{\\partial y^2} = 0 ∂x2∂2p​+∂y2∂2p​=0
2. 通过观察这个方程可以发现，两项是关于x和y的扩散项方程，因此可以用二阶中心差分进行离散：
   p i + 1 , j n − 2 p i , j n + p i − 1 , j n Δ x 2 + p i , j + 1 n − 2 p i , j n + p i , j − 1 n Δ y 2 = 0 \\frac{p_{i+1, j}^n - 2p_{i,j}^n + p_{i-1,j}^n}{\\Delta x^2} + \\frac{p_{i,j+1}^n - 2p_{i,j}^n + p_{i, j-1}^n}{\\Delta y^2} = 0 Δx2pi+1,jn​−2pi,jn​+pi−1,jn​​+Δy2pi,j+1n​−2pi,jn​+pi,j−1n​​=0
3. 该方程与时间t无关（也就是常说的稳态），那么可以通过迭代的方法求解 p i , j n p_{i,j}^n pi,jn​ ，即将离散方程转化为五点差分形式：
   p i , j n = Δ y 2 ( p i + 1 , j n + p i − 1 , j n ) + Δ x 2 ( p i , j + 1 n + p i , j − 1 n ) 2 ( Δ x 2 + Δ y 2 ) p_{i,j}^n = \\frac{\\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\\Delta x^2 + \\Delta y^2)} pi,jn​=2(Δx2+Δy2)Δy2(pi+1,jn​+pi−1,jn​)+Δx2(pi,j+1n​+pi,j−1n​)​
4. 初始条件：$p=0$
5. 边界条件：
   • $p=0$ at $x=0$
   • $p=y$ at $x=2$
   • ∂ p ∂ y = 0 \\frac{\\partial p}{\\partial y}=0 ∂y∂p​=0 at $y=0,1$
6. 解析解：
   p ( x , y ) = x 4 − 4 ∑ n = 1 , o d d ∞ 1 ( n π ) 2 sinh ⁡ 2 n π sinh ⁡ n π x cos ⁡ n π y p(x,y)=\\frac{x}{4}-4\\sum_{n=1,odd}^{\\infty}\\frac{1}{(n\\pi)^2\\sinh2n\\pi}\\sinh n\\pi x\\cos n\\pi y p(x,y)=4x​−4n=1,odd∑∞​(nπ)2sinh2nπ1​sinhnπxcosnπy
7. 代码如下（自己写的，原文可以看最上方的链接）

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cmdef plot2D(x,y,p):plt.ion() # 这个东西是方便close，不然close不鸟fig=plt.figure(figsize=(11,7),dpi=100) # 在Matplotlib 3.4版本之后，将fig.gca弃用了ax = fig.add_subplot(projection='3d')X,Y=np.meshgrid(x,y)surf=ax.plot_surface(X,Y,p[:],rstride=1,cstride=1,cmap=cm.viridis,linewidth=0,antialiased=False)ax.set_xlim(0,2)ax.set_ylim(0, 1)ax.view_init(30,225)# create x and y coordinate
nx=100
ny=100
x=np.linspace(0,2,nx)
y=np.linspace(0,1,ny)
dx=2/(nx-1)
dy=1/(ny-1)# initial p
p=np.zeros((nx,ny))# boundary p
p[0,:]=0
p[-1,:]=y
p[:,0]=p[:,1] # dp/dy=0 at y=0
p[:,-1]=p[:,-2] # dp/dy=0 at y=1# start iteration
error=1
pn=np.empty_like(p)
while error>=1e-6:pn=p.copy()p[1:-1,1:-1]=((dy2*(p[2:,1:-1]+p[0:-2,1:-1]))+dx2*(p[1:-1,2:]+p[1:-1,0:-2]))/(2*(dx2+dy2))# boundaryp[0, :] = 0p[-1, :] = yp[:, 0] = p[:, 1]  # dp/dy=0 at y=0p[:, -1] = p[:, -2]  # dp/dy=0 at y=1# error= norm L2error=np.linalg.norm(x=p-pn,ord=2)# draw in every iterationplt.close()plot2D(x,y,p)plt.show()plt.pause(1)

Step 10: 2D Poisson Equation

1. 这一步作者在最后提到，其实求解这种方程的代码都很像（事实也确实如此，仅仅是边界条件和初始条件不同，离散方程的方法大差不差），如果想规整这些代码，涉及到python中package的概念。
2. 二维泊松方程的形式如下：
   ∂ 2 p ∂ x 2 + ∂ 2 p ∂ y 2 = b \\frac{\\partial ^2 p}{\\partial x^2} + \\frac{\\partial ^2 p}{\\partial y^2} = b ∂x2∂2p​+∂y2∂2p​=b
3. 离散方式和上文大同小异，就是扩散项和常数源项：
   p i + 1 , j n − 2 p i , j n + p i − 1 , j n Δ x 2 + p i , j + 1 n − 2 p i , j n + p i , j − 1 n Δ y 2 = b i , j n \\frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\\Delta x^2}+\\frac{p_{i,j+1}^{n}-2 p_{i,j}^{n}+p_{i,j-1}^{n}}{\\Delta y^2}=b_{i,j}^{n} Δx2pi+1,jn​−2pi,jn​+pi−1,jn​​+Δy2pi,j+1n​−2pi,jn​+pi,j−1n​​=bi,jn​
4. 代码就不贴了，有需要的话看原文去，这里也让笔者大受启发，在求解CFD问题的时候，不要想着一口气吃成一个大胖子，先从稳态方程写起，弄好之后再加非稳态项，最后加源项。

Step 11: Cavity Flow with Navier–Stokes

1. 这是N-S方程质量和动量守恒方程组：
   ∇ ⋅ v ⃗ = 0 \\nabla\\cdot\\vec{v} = 0 ∇⋅v =0
   ∂ v ⃗ ∂ t + ( v ⃗ ⋅ ∇ ) v ⃗ = − 1 ρ ∇ p + ν ∇ 2 v ⃗ \\frac{\\partial \\vec{v}}{\\partial t}+(\\vec{v}\\cdot\\nabla)\\vec{v} = -\\frac{1}{\\rho}\\nabla p + \\nu \\nabla^2\\vec{v} ∂t∂v ​+(v ⋅∇)v =−ρ1​∇p+ν∇2v
2. 二维方腔流方程形式如下，其中第三个压力的扩散方程就是上一步推出的2维泊松方程，这里由于要分别求压力和流速，可以尝试采用SIMPLE算法
   ∂ u ∂ t + u ∂ u ∂ x + v ∂ u ∂ y = − 1 ρ ∂ p ∂ x + ν ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) \\frac{\\partial u}{\\partial t}+u\\frac{\\partial u}{\\partial x}+v\\frac{\\partial u}{\\partial y} = -\\frac{1}{\\rho}\\frac{\\partial p}{\\partial x}+\\nu \\left(\\frac{\\partial^2 u}{\\partial x^2}+\\frac{\\partial^2 u}{\\partial y^2} \\right) ∂t∂u​+u∂x∂u​+v∂y∂u​=−ρ1​∂x∂p​+ν(∂x2∂2u​+∂y2∂2u​)

∂ v ∂ t + u ∂ v ∂ x + v ∂ v ∂ y = − 1 ρ ∂ p ∂ y + ν ( ∂ 2 v ∂ x 2 + ∂ 2 v ∂ y 2 ) \\frac{\\partial v}{\\partial t}+u\\frac{\\partial v}{\\partial x}+v\\frac{\\partial v}{\\partial y} = -\\frac{1}{\\rho}\\frac{\\partial p}{\\partial y}+\\nu\\left(\\frac{\\partial^2 v}{\\partial x^2}+\\frac{\\partial^2 v}{\\partial y^2}\\right) ∂t∂v​+u∂x∂v​+v∂y∂v​=−ρ1​∂y∂p​+ν(∂x2∂2v​+∂y2∂2v​)

∂ 2 p ∂ x 2 + ∂ 2 p ∂ y 2 = − ρ ( ∂ u ∂ x ∂ u ∂ x + 2 ∂ u ∂ y ∂ v ∂ x + ∂ v ∂ y ∂ v ∂ y ) \\frac{\\partial^2 p}{\\partial x^2}+\\frac{\\partial^2 p}{\\partial y^2} = -\\rho\\left(\\frac{\\partial u}{\\partial x}\\frac{\\partial u}{\\partial x}+2\\frac{\\partial u}{\\partial y}\\frac{\\partial v}{\\partial x}+\\frac{\\partial v}{\\partial y}\\frac{\\partial v}{\\partial y} \\right) ∂x2∂2p​+∂y2∂2p​=−ρ(∂x∂u​∂x∂u​+2∂y∂u​∂x∂v​+∂y∂v​∂y∂v​)

1. 方程离散，作者建议这里要手写一下，原理不难，就是简单的差分罢了，但是很考验一个人的手感：
   u i , j n + 1 − u i , j n Δ t + u i , j n u i , j n − u i − 1 , j n Δ x + v i , j n u i , j n − u i , j − 1 n Δ y = − 1 ρ p i + 1 , j n − p i − 1 , j n 2 Δ x + ν ( u i + 1 , j n − 2 u i , j n + u i − 1 , j n Δ x 2 + u i , j + 1 n − 2 u i , j n + u i , j − 1 n Δ y 2 ) \\begin{split} & \\frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\\Delta t}+u_{i,j}^{n}\\frac{u_{i,j}^{n}-u_{i-1,j}^{n}}{\\Delta x}+v_{i,j}^{n}\\frac{u_{i,j}^{n}-u_{i,j-1}^{n}}{\\Delta y} = \\\\ & \\qquad -\\frac{1}{\\rho}\\frac{p_{i+1,j}^{n}-p_{i-1,j}^{n}}{2\\Delta x}+\\nu\\left(\\frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{\\Delta x^2}+\\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{\\Delta y^2}\\right) \\end{split} ​Δtui,jn+1​−ui,jn​​+ui,jn​Δxui,jn​−ui−1,jn​​+vi,jn​Δyui,jn​−ui,j−1n​​=−ρ1​2Δxpi+1,jn​−pi−1,jn​​+ν(Δx2ui+1,jn​−2ui,jn​+ui−1,jn​​+Δy2ui,j+1n​−2ui,jn​+ui,j−1n​​)​
   v i , j n + 1 − v i , j n Δ t + u i , j n v i , j n − v i − 1 , j n Δ x + v i , j n v i , j n − v i , j − 1 n Δ y = − 1 ρ p i , j + 1 n − p i , j − 1 n 2 Δ y + ν ( v i + 1 , j n − 2 v i , j n + v i − 1 , j n Δ x 2 + v i , j + 1 n − 2 v i , j n + v i , j − 1 n Δ y 2 ) \\begin{split} &\\frac{v_{i,j}^{n+1}-v_{i,j}^{n}}{\\Delta t}+u_{i,j}^{n}\\frac{v_{i,j}^{n}-v_{i-1,j}^{n}}{\\Delta x}+v_{i,j}^{n}\\frac{v_{i,j}^{n}-v_{i,j-1}^{n}}{\\Delta y} = \\\\ & \\qquad -\\frac{1}{\\rho}\\frac{p_{i,j+1}^{n}-p_{i,j-1}^{n}}{2\\Delta y} +\\nu\\left(\\frac{v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}}{\\Delta x^2}+\\frac{v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}}{\\Delta y^2}\\right) \\end{split} ​Δtvi,jn+1​−vi,jn​​+ui,jn​Δxvi,jn​−vi−1,jn​​+vi,jn​Δyvi,jn​−vi,j−1n​​=−ρ1​2Δypi,j+1n​−pi,j−1n​​+ν(Δx2vi+1,jn​−2vi,jn​+vi−1,jn​​+Δy2vi,j+1n​−2vi,jn​+vi,j−1n​​)​
   p i + 1 , j n − 2 p i , j n + p i − 1 , j n Δ x 2 + p i , j + 1 n − 2 p i , j n + p i , j − 1 n Δ y 2 = ρ [ 1 Δ t ( u i + 1 , j − u i − 1 , j 2 Δ x + v i , j + 1 − v i , j − 1 2 Δ y ) − u i + 1 , j − u i − 1 , j 2 Δ x u i + 1 , j − u i − 1 , j 2 Δ x − 2 u i , j + 1 − u i , j − 1 2 Δ y v i + 1 , j − v i − 1 , j 2 Δ x − v i , j + 1 − v i , j − 1 2 Δ y v i , j + 1 − v i , j − 1 2 Δ y ] \\begin{split} & \\frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\\Delta x^2}+\\frac{p_{i,j+1}^{n}-2p_{i,j}^{n}+p_{i,j-1}^{n}}{\\Delta y^2} = \\\\ & \\qquad \\rho \\left[ \\frac{1}{\\Delta t}\\left(\\frac{u_{i+1,j}-u_{i-1,j}}{2\\Delta x}+\\frac{v_{i,j+1}-v_{i,j-1}}{2\\Delta y}\\right) -\\frac{u_{i+1,j}-u_{i-1,j}}{2\\Delta x}\\frac{u_{i+1,j}-u_{i-1,j}}{2\\Delta x} - 2\\frac{u_{i,j+1}-u_{i,j-1}}{2\\Delta y}\\frac{v_{i+1,j}-v_{i-1,j}}{2\\Delta x} - \\frac{v_{i,j+1}-v_{i,j-1}}{2\\Delta y}\\frac{v_{i,j+1}-v_{i,j-1}}{2\\Delta y}\\right] \\end{split} ​Δx2pi+1,jn​−2pi,jn​+pi−1,jn​​+Δy2pi,j+1n​−2pi,jn​+pi,j−1n​​=ρ[Δt1​(2Δxui+1,j​−ui−1,j​​+2Δyvi,j+1​−vi,j−1​​)−2Δxui+1,j​−ui−1,j​​2Δxui+1,j​−ui−1,j​​−22Δyui,j+1​−ui,j−1​​2Δxvi+1,j​−vi−1,j​​−2Δyvi,j+1​−vi,j−1​​2Δyvi,j+1​−vi,j−1​​]​
2. 把上一时刻的压力和速度放一侧：
   u i , j n + 1 = u i , j n − u i , j n Δ t Δ x ( u i , j n − u i − 1 , j n ) − v i , j n Δ t Δ y ( u i , j n − u i , j − 1 n ) − Δ t ρ 2 Δ x ( p i + 1 , j n − p i − 1 , j n ) + ν ( Δ t Δ x 2 ( u i + 1 , j n − 2 u i , j n + u i − 1 , j n ) + Δ t Δ y 2 ( u i , j + 1 n − 2 u i , j n + u i , j − 1 n ) ) \\begin{split} u_{i,j}^{n+1} = u_{i,j}^{n} & - u_{i,j}^{n} \\frac{\\Delta t}{\\Delta x} \\left(u_{i,j}^{n}-u_{i-1,j}^{n}\\right) - v_{i,j}^{n} \\frac{\\Delta t}{\\Delta y} \\left(u_{i,j}^{n}-u_{i,j-1}^{n}\\right) \\\\ & - \\frac{\\Delta t}{\\rho 2\\Delta x} \\left(p_{i+1,j}^{n}-p_{i-1,j}^{n}\\right) \\\\ & + \\nu \\left(\\frac{\\Delta t}{\\Delta x^2} \\left(u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}\\right) + \\frac{\\Delta t}{\\Delta y^2} \\left(u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}\\right)\\right) \\end{split} ui,jn+1​=ui,jn​​−ui,jn​ΔxΔt​(ui,jn​−ui−1,jn​)−vi,jn​ΔyΔt​(ui,jn​−ui,j−1n​)−ρ2ΔxΔt​(pi+1,jn​−pi−1,jn​)+ν(Δx2Δt​(ui+1,jn​−2ui,jn​+ui−1,jn​)+Δy2Δt​(ui,j+1n​−2ui,jn​+ui,j−1n​))​
   p i , j n = ( p i + 1 , j n + p i − 1 , j n ) Δ y 2 + ( p i , j + 1 n + p i , j − 1 n ) Δ x 2 2 ( Δ x 2 + Δ y 2 ) − ρ Δ x 2 Δ y 2 2 ( Δ x 2 + Δ y 2 ) × [ 1 Δ t ( u i + 1 , j − u i − 1 , j 2 Δ x + v i , j + 1 − v i , j − 1 2 Δ y ) − u i + 1 , j − u i − 1 , j 2 Δ x u i + 1 , j − u i − 1 , j 2 Δ x − 2 u i , j + 1 − u i , j − 1 2 Δ y v i + 1 , j − v i − 1 , j 2 Δ x − v i , j + 1 − v i , j − 1 2 Δ y v i , j + 1 − v i , j − 1 2 Δ y ] \\begin{split} p_{i,j}^{n} = & \\frac{\\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\\right) \\Delta y^2 + \\left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\\right) \\Delta x^2}{2\\left(\\Delta x^2+\\Delta y^2\\right)} \\\\ & -\\frac{\\rho\\Delta x^2\\Delta y^2}{2\\left(\\Delta x^2+\\Delta y^2\\right)} \\\\ & \\times \\left[\\frac{1}{\\Delta t}\\left(\\frac{u_{i+1,j}-u_{i-1,j}}{2\\Delta x}+\\frac{v_{i,j+1}-v_{i,j-1}}{2\\Delta y}\\right)-\\frac{u_{i+1,j}-u_{i-1,j}}{2\\Delta x}\\frac{u_{i+1,j}-u_{i-1,j}}{2\\Delta x} -2\\frac{u_{i,j+1}-u_{i,j-1}}{2\\Delta y}\\frac{v_{i+1,j}-v_{i-1,j}}{2\\Delta x}-\\frac{v_{i,j+1}-v_{i,j-1}}{2\\Delta y}\\frac{v_{i,j+1}-v_{i,j-1}}{2\\Delta y}\\right] \\end{split} pi,jn​=​2(Δx2+Δy2)(pi+1,jn​+pi−1,jn​)Δy2+(pi,j+1n​+pi,j−1n​)Δx2​−2(Δx2+Δy2)ρΔx2Δy2​×[Δt1​(2Δxui+1,j​−ui−1,j​​+2Δyvi,j+1​−vi,j−1​​)−2Δxui+1,j​−ui−1,j​​2Δxui+1,j​−ui−1,j​​−22Δyui,j+1​−ui,j−1​​2Δxvi+1,j​−vi−1,j​​−2Δyvi,j+1​−vi,j−1​​2Δyvi,j+1​−vi,j−1​​]​
3. 边界条件：

   • $u=1$ at $y=2$ (the “lid”);

   • $u,v=0$ on the other boundaries;

   • ∂ p ∂ y = 0 \\frac{\\partial p}{\\partial y}=0 ∂y∂p​=0 at $y=0$;

   • $p=0$ at $y=2$

   • ∂ p ∂ x = 0 \\frac{\\partial p}{\\partial x}=0 ∂x∂p​=0 at $x=0,2$

4. 代码如下：

import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inlinenx = 41
ny = 41
nt = 500
nit = 50
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 2, ny)
X, Y = numpy.meshgrid(x, y)rho = 1
nu = .1
dt = .001u = numpy.zeros((ny, nx))
v = numpy.zeros((ny, nx))
p = numpy.zeros((ny, nx)) 
b = numpy.zeros((ny, nx))def build_up_b(b, rho, dt, u, v, dx, dy):b[1:-1, 1:-1] = (rho * (1 / dt * ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) + (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))2 -2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *(v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx))-((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))2))return bdef pressure_poisson(p, dx, dy, b):pn = numpy.empty_like(p)pn = p.copy()for q in range(nit):pn = p.copy()p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy2 + (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx2) /(2 * (dx2 + dy2)) -dx2 * dy2 / (2 * (dx2 + dy2)) * b[1:-1,1:-1])p[:, -1] = p[:, -2] # dp/dx = 0 at x = 2p[0, :] = p[1, :]   # dp/dy = 0 at y = 0p[:, 0] = p[:, 1]   # dp/dx = 0 at x = 0p[-1, :] = 0        # p = 0 at y = 2return pdef cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu):un = numpy.empty_like(u)vn = numpy.empty_like(v)b = numpy.zeros((ny, nx))for n in range(nt):un = u.copy()vn = v.copy()b = build_up_b(b, rho, dt, u, v, dx, dy)p = pressure_poisson(p, dx, dy, b)u[1:-1, 1:-1] = (un[1:-1, 1:-1]-un[1:-1, 1:-1] * dt / dx *(un[1:-1, 1:-1] - un[1:-1, 0:-2]) -vn[1:-1, 1:-1] * dt / dy *(un[1:-1, 1:-1] - un[0:-2, 1:-1]) -dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, 0:-2]) +nu * (dt / dx2 *(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +dt / dy2 *(un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])))v[1:-1,1:-1] = (vn[1:-1, 1:-1] -un[1:-1, 1:-1] * dt / dx *(vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -vn[1:-1, 1:-1] * dt / dy *(vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -dt / (2 * rho * dy) * (p[2:, 1:-1] - p[0:-2, 1:-1]) +nu * (dt / dx2 *(vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +dt / dy2 *(vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))u[0, :]  = 0u[:, 0]  = 0u[:, -1] = 0u[-1, :] = 1    # set velocity on cavity lid equal to 1v[0, :]  = 0v[-1, :] = 0v[:, 0]  = 0v[:, -1] = 0return u, v, pu = numpy.zeros((ny, nx))
v = numpy.zeros((ny, nx))
p = numpy.zeros((ny, nx))
b = numpy.zeros((ny, nx))
nt = 100 # change to 700
u, v, p = cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)fig = pyplot.figure(figsize=(11,7), dpi=100)
# plotting the pressure field as a contour
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)  
pyplot.colorbar()
# plotting the pressure field outlines
pyplot.contour(X, Y, p, cmap=cm.viridis)  
# plotting velocity field
pyplot.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2]) 
pyplot.xlabel('X')
pyplot.ylabel('Y');

4. 也可以用streamplot画图，这样流线是连续的，用quiver画图流线是间断的。

fig = pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
pyplot.colorbar()
pyplot.contour(X, Y, p, cmap=cm.viridis)
pyplot.streamplot(X, Y, u, v)
pyplot.xlabel('X')
pyplot.ylabel('Y');

Step 12: Channel Flow with Navier–Stokes

1. 二维管道流就是在u方向的动量方程加了源项F，原文处是采用迭代收敛的方法计算的
2. 全篇看完后第一感觉就是太浅显了，全文基本是通过有限差分的方式进行求解，然后对于其物理本质基本没有概括。但总而言之，浅显有浅显的好处，非常利好初学者自学，但笔者后期时间有限，加上作者这最后几个步骤同质化过于严重，导致笔者第三篇笔记写的非常水，所以强烈建议有兴趣的读者去看看github原文