基于Python的CFD编程入门（3）一维扩散方程
讨论过了基础的对流问题，接下来介绍一下扩散（Diffusion）问题。
对流和扩散的区别在于：
对流是流体微团宏观的定向运动，带有强烈的方向性。在对流的作用下， 发生在某一地点上的扰动只能向其下游方向传递而不会逆向传播。
扩散是由分子的不规则热运动所致，分子不规则热运动对空间不同方向的几率都是一样的，因而扩散过程可以把发生在某一地点上的扰动的影响向各个方向传递。
理论基础
一维扩散方程为：
由于扩散方程为二阶偏导，可以使用中心差分的方法进行离散化。
中心差分由一阶导的前向差分和后向差分之和得到:
通过中心差分法离散扩散方程：
则需要求解的未知数：
代码实现
这里依旧使用和求解对流问题时同样的初始化方法：
import numpy as npimport matplotlib.pyplot as pltnx = 41dx = 2 / (nx - 1)nt = 20nu = 0.3sigma = .2 #这里的参数后续会介绍到dt = sigma * dx**2 / nuu = np.ones(nx)u[int(.5 / dx) : int(1 / dx + 1)] = 2un = np.ones(nx)grid = np.linspace(0,2,nx)plt.plot(grid,u)
迭代部分和对流方程类似，注意空间点的取值范围即可：
for n in range(nt):un = u.copy()for i in range(1, nx-1):u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2*un[i] + un[i-1])plt.plot(grid,u)
计算结果为：
每个时间步的结果为：
方程的稳定性
既然用到了中心差分的方法，那么我们就不得不考虑计算的稳定性。
当我们尝试增加粘性、空间点数或者时间步长时，在一定情况下中心差分法就会变得不稳定。
这是由比率r决定，如果 r 太大，该方法就会变得不稳定：
1.初始值
2.增加粘度
3.增加空间点
4.增加计算步长
关于稳定性后续会进行专题讨论，本期只做简单介绍。
计算中常见的三种不稳定性：
初值问题显式格式的不稳定性：显式格式求解抛物型方程，由于时间步长取得过大，引起解的振荡，失去物理意义。这时需要进一步确定计算时间步长。
代数方程求解过程中出现的不稳定性：迭代收敛条件不满足，导致迭代求解不收敛，这是求解过程的不稳定性问题（得不到解）。
对流项离散格式引起的不稳定性：采用CD（中心差分），TUD（三阶迎风），QUICK（对流项的二次迎风插值）求解稳态对流问题时，如流速很高，网格较粗也会产生振荡的解，称为对流项离散格式的不稳态定性。此时需要找出产生振荡的临界Peclet数
对于大多数现有的对流-扩散方程的离散格式一般采用这样的方式来研究其稳定性：即把所分析的离散格式应用于一维稳态无源项的模型方程，通过对这一离散方程稳定性的分析，找出该格式的稳定性条件，即找出使数值解不出现振荡的网格Peclet（佩克莱）数的极限值（称为临界Peclet数）。
正系数法
离散方程精确解分析法
反馈灵敏度分析法
“符号不变”原则（sign preservation rule）
