基于Python的CFD编程入门（1）— 一维线性对流问题
我们知道求解NS方程是CFD的核心问题，即便平时使用的是商业软件，通过学习和编写代码求解CFD相关方程可以帮助我们了解CFD仿真软件的底层逻辑与误差控制，让我们对CFD的理解更上一层楼。
虽然Python的运行效率不足以让我们在此基础上开发CFD求解器，但是其简单易用的特点可以让我们把精力聚焦于问题本身，能够更快的对CFD算法等问题有清晰的认识。
理论基础
不可压动量方程：
我们只保留速度的x方向分量，再令速度的一阶导系数为常数，这样就得到了一维线性对流（1D Linear Convection）方程。
该方程表示波的传播没有形状变化，而速度为常系数c。
给定 t = 0 时初始条件 u(x,0)，则解为u(x,0)。
所以方程的解为：
可以认为是波以匀速c沿着x轴正方向平移。
编写程序求解上述偏微分方程，就需要在空间和时间上离散这个方程，这里对于时间步长，使用上标 n，对于空间步长（网格单元），使用下标 i，所以对时间导数采用前向差分有：
对空间导数采用后向差分方案：
则一维对流方程的离散形式为：
代码实现
使用Python实现上述方程的求解，这里需要使用到numpy、matplotlib库。
由于是一维问题，所以其计算域为一条直线，我们令计算域长度为2。
初始化时，令所有位置的速度都是0，而（0.5，1）的位置速度为1。
令直线上的网格单元数为nx ，单个尺寸为dx ，时间步数量为nt ，时间步长为dt 。
通过下述代码表示初始时刻的速度分布
import numpy as npimport matplotlib.pyplot as pltimport time,sysnx = 41 # 网格数dx = 2 / (nx-1)nt = 30 #时间步数dt = 0.025 #时间步长c = 1u = np.zeros(nx)u[int(.5/dx):int(1/dx + 1)] = 1plt.plot(np.linspace(0, 2, nx),u)
图像表示为：
上图类似方波的图像也可以称作"帽函数(hat function)",其表示了各个点初始时刻的速度u.
可以注意到两条竖边并没有完全垂直，这是由网格单元的尺寸造成的。当我们划分的数量越多，这两条边就越接近竖直。
接下来需要编写两个循环来分别实现位置和时间的迭代：
un = np.ones(nx)grid = np.linspace(0, 2, nx)for n in range(nt):un = u.copy()for i in range(1, nx):u[i] = un[i] - c * dt / dx * (un[i] - un[i - 1])plt.plot(grid, u)
迭代计算后的结果为：
这里打印出每个时间步的计算结果可以更加直观的理解:
很明显，波形向正方向传播，形状没有变化；波形的最高点逐步降低，直到解达到收敛。
增加计算域长度和时间步数量，可以发现波形的振幅有衰减的趋势，也就是速度会逐渐衰减。
