最近感觉CSDN有点大病,用各种浏览器都不显示最上面那一栏(就是点赞、数据、发布啥啥的那一栏),今天用ubuntu偶然发现又显示了, 赶紧把之前写的东西发出来记录一下,不知道问题出在哪:(
CFDpython - 12 steps to N-S equation
后面四步,每一步都是一个不同的方程,分别是:二维拉普拉斯方程、二维泊松方程、二维空腔流动、二维管渠流动
Step 9: 2D Laplace Equation
- 方程形式如下:
∂ 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 - 通过观察这个方程可以发现,两项是关于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 - 该方程与时间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) - 初始条件: p = 0 p=0 p=0
- 边界条件:
- p = 0 p=0 p=0 at x = 0 x=0 x=0
- p = y p=y p=y at x = 2 x=2 x=2
- ∂ p ∂ y = 0 \frac{\partial p}{\partial y}=0 ∂y∂p=0 at y = 0 , 1 y=0, \ 1 y=0, 1
- 解析解:
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π1sinhnπxcosnπy - 代码如下(自己写的,原文可以看最上方的链接)
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
def 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]=((dy**2*(p[2:,1:-1]+p[0:-2,1:-1]))+dx**2*(p[1:-1,2:]+p[1:-1,0:-2]))/(2*(dx**2+dy**2))
# boundary
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
# error= norm L2
error=np.linalg.norm(x=p-pn,ord=2)
# draw in every iteration
plt.close()
plot2D(x,y,p)
plt.show()
plt.pause(1)
Step 10: 2D Poisson Equation
- 这一步作者在最后提到,其实求解这种方程的代码都很像(事实也确实如此,仅仅是边界条件和初始条件不同,离散方程的方法大差不差),如果想规整这些代码,涉及到python中
package
的概念。 - 二维泊松方程的形式如下:
∂ 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 - 离散方式和上文大同小异,就是扩散项和常数源项:
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 - 代码就不贴了,有需要的话看原文去,这里也让笔者大受启发,在求解CFD问题的时候,不要想着一口气吃成一个大胖子,先从稳态方程写起,弄好之后再加非稳态项,最后加源项。
Step 11: Cavity Flow with Navier–Stokes
- 这是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维泊松方程,这里由于要分别求压力和流速,可以尝试采用
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)
- 方程离散,作者建议这里要手写一下,原理不难,就是简单的差分罢了,但是很考验一个人的手感:
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=−ρ12Δ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=−ρ12Δ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,j2Δxui+1,j−ui−1,j−22Δyui,j+1−ui,j−12Δxvi+1,j−vi−1,j−2Δyvi,j+1−vi,j−12Δyvi,j+1−vi,j−1] - 把上一时刻的压力和速度放一侧:
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,j2Δxui+1,j−ui−1,j−22Δyui,j+1−ui,j−12Δxvi+1,j−vi−1,j−2Δyvi,j+1−vi,j−12Δyvi,j+1−vi,j−1] - 边界条件:
-
u = 1 u=1 u=1 at y = 2 y=2 y=2 (the “lid”);
-
u , v = 0 u, v=0 u,v=0 on the other boundaries;
-
∂ p ∂ y = 0 \frac{\partial p}{\partial y}=0 ∂y∂p=0 at y = 0 y=0 y=0;
-
p = 0 p=0 p=0 at y = 2 y=2 y=2
-
∂ p ∂ x = 0 \frac{\partial p}{\partial x}=0 ∂x∂p=0 at x = 0 , 2 x=0,2 x=0,2
-
- 代码如下:
import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
nx = 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 = .001
u = 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 b
def 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]) * dy**2 +
(pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /
(2 * (dx**2 + dy**2)) -
dx**2 * dy**2 / (2 * (dx**2 + dy**2)) *
b[1:-1,1:-1])
p[:, -1] = p[:, -2] # dp/dx = 0 at x = 2
p[0, :] = p[1, :] # dp/dy = 0 at y = 0
p[:, 0] = p[:, 1] # dp/dx = 0 at x = 0
p[-1, :] = 0 # p = 0 at y = 2
return p
def 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 / dx**2 *
(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
dt / dy**2 *
(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 / dx**2 *
(vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
dt / dy**2 *
(vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))
u[0, :] = 0
u[:, 0] = 0
u[:, -1] = 0
u[-1, :] = 1 # set velocity on cavity lid equal to 1
v[0, :] = 0
v[-1, :] = 0
v[:, 0] = 0
v[:, -1] = 0
return u, v, p
u = 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');
- 也可以用
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
- 二维管道流就是在u方向的动量方程加了源项F,原文处是采用迭代收敛的方法计算的
- 全篇看完后第一感觉就是太浅显了,全文基本是通过有限差分的方式进行求解,然后对于其物理本质基本没有概括。但总而言之,浅显有浅显的好处,非常利好初学者自学,但笔者后期时间有限,加上作者这最后几个步骤同质化过于严重,导致笔者第三篇笔记写的非常水,所以强烈建议有兴趣的读者去看看github原文