【Python4CFD】笔记step9-12

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

CFDpython - 12 steps to N-S equation

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

Step 9: 2D Laplace Equation

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

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

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

Step 11: Cavity Flow with Navier–Stokes

1. 这是N-S方程质量和动量守恒方程组：
   $$\nabla \cdot \vec{v}=0$$
   $$\frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot \nabla )\vec{v}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\vec{v}$$
2. 二维方腔流方程形式如下，其中第三个压力的扩散方程就是上一步推出的2维泊松方程，这里由于要分别求压力和流速，可以尝试采用SIMPLE算法
   $$\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)$$

$$\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)$$

$$\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)$$

1. 方程离散，作者建议这里要手写一下，原理不难，就是简单的差分罢了，但是很考验一个人的手感：
   $$\begin{array}{rl}& \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}=\\ & -\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{array}$$
   $$\begin{array}{rl}& \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}=\\ & -\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{array}$$
   $$\begin{array}{rl}& \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}=\\ & \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{array}$$
2. 把上一时刻的压力和速度放一侧：
   $$\begin{array}{rl}{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{array}$$
   $$\begin{array}{rl}{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{array}$$
3. 边界条件：

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

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

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

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

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

4. 代码如下：

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');

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原文