2. 扩散方程

文档B站视频讲解：

LBM入门课程(2.1)LBM求解扩散方程(一维)

LBM入门课程(2.2)LBM的时空二阶精度及变系数/带源项的扩散方程求解

LBM入门课程(2.3)D2Q5/D2Q9求扩散方程及非平衡态外推边界条件

这节我们讲解扩散方程的格子玻尔兹曼方法的实现。包括一维和二维的情况。

2.1. 什么是扩散方程？

扩散方程是描述物质扩散的方程。在一维情况下，扩散方程是这样的：

\[ \frac{\partial \phi}{\partial t} = D \frac{\partial^2 \phi}{\partial x^2} \]

其中\(\phi\)可以是温度/物质浓度/速度等，\(D\)是扩散系数。在高维情况下，扩散方程是这样的：

\[ \frac{\partial \phi}{\partial t} = D \nabla^2 \phi \]

很多人接触流体力学和LBM时，之所以头大，就是因为看不懂上面的这种符号，我这里给大家展开一下。\(\nabla^2\)是拉普拉斯算子，它在直角坐标系下的表达式是：

\[ \nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} \]

因此，三维情况下的扩散方程就是：

\[ \frac{\partial \phi}{\partial t} = D \left( \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} + \frac{\partial^2 \phi}{\partial z^2} \right) \]

扩散方程的直观算例，下面这个图是左边固定一个\(T=1\)的边界条件，右边固定一个\(T=0\)的边界条件，然后就会慢慢的向右边扩散。这就是称它为扩散方程的原因。

下图则展示了一个初始条件为正弦分布场随时间的演化，两侧为周期性边界条件。

相信大家通过这两个算例已经对扩散方程有了一个直观的认识。接下来我们就来看看如何用LBM来求解扩散方程的。

2.2. 一维扩散方程的D1Q2模型

重新写一个一维扩散方程：

\[ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} \]

当把\(T\)当作变量时，\(\alpha\)一般称为热扩散系数。

2.2.1. D1Q2模型整理：

两个离散速度

\[ e_1 = c, \quad e_2 = -c \]

其中\(c=\delta_x/\delta_t\), \(\delta_x\)是空间步长，\(\delta_t\)是时间步长。

对应的权重

\[ w_1 = w_2 = 0.5 \]

平衡态分布函数

\[ f_i^\mathrm{eq} = w_i T \]

演化方程

\[ f_i(\boldsymbol{x}+\boldsymbol{e}_i \delta_t,t+\delta_t)=f_i(\boldsymbol{x},t)-\frac{1}{\tau}[f_i(\boldsymbol{x},t)-f_i^{\mathrm{eq}}(\boldsymbol{x},t)] \]

也可以写成如下的形式：

\[ f_i(\boldsymbol{x}+\boldsymbol{e}_i \delta_t,t+\delta_t)=(1-\omega) f_i(\boldsymbol{x},t)+\omega f_i^{\mathrm{eq}}(\boldsymbol{x},t) \]

其中\(\omega = 1/\tau\)。

宏观量

\[ T= \sum_i f_i \]

热扩散率

\[ \alpha = \frac{\delta _x^2}{\delta_t}\left(\tau - \frac{1}{2}\right) \]

2.2.2. 边界条件

边界条件是求解扩散方程的关键。在一维情况下，我们一般会有两种边界条件：

2.2.2.1. 第一类边界条件：恒温边界条件

这种边界条件是指在某个位置，温度是恒定的。比如在\(x=0\)处，\(T=1\)，在\(x=L\)处，\(T=0\)。这种边界条件的处理方法是：

左壁面

\[ f_1 = T_\mathrm{left} - f_2 \]

右壁面

\[ f_2 = T_\mathrm{right} - f_1 \]

2.2.2.2. 第二类边界条件：恒热流边界条件

这种边界条件是指在某个位置，热流密度是恒定的。比如在\(x=0\)处，$\(\frac{\partial T}{\partial x} = q\)$ 对于这种边界条件，我们需要知道的是，在D1Q2中，每个地方的热流密度可以通过以下方式计算：

\[ q = k\omega\frac{f_1 - f_2}{\delta_x} \]

其中，\(k\)为热导率。因此，我们可以通过以下方式来处理边界条件：

左壁面

\[ f_1 = f_2 + \frac{q\delta_x}{k\omega} \]

右壁面

\[ f_2 = f_1 - \frac{q\delta_x}{k\omega} \]

2.2.3. 算例验证

2.2.3.1. 恒温边界条件

左侧为\(T=1\)，右侧为\(T=0\)。计算域长为\(100\)，初始条件为\(T=0.5\)。计算的参数设置为热扩散率\(\alpha=0.25\)，\(\delta_x=1\)，\(\delta_t=1\)。我们可以看到，随着时间的推移，温度会慢慢的变成线性。

绘制热流密度的分布，FDM表示有限差分，计算方法为：

\[ q = -k \frac{\partial T}{\partial x} = -k \frac{T_{i+1}-T_{i}}{\delta_x} \]

LBM的计算方法为：

\[ q = k\omega \frac{f_1 - f_2}{\delta_x} \]

计算时，我们取\(k=1\)。可以看到两种方式计算的结果是一致的。

下面展示的是\(\delta_x = 0.5\)，\(\delta_t=0.5\)时的计算结果：

可以看到我们的程序可以独立的调节\(\delta_x\)和\(\delta_t\)，并且计算结果是正确的。

2.2.3.2. 恒热流边界条件

此算例的边界条件为左侧恒定热流边界条件：\(q=0.01\)，右侧为恒定壁温边界条件：\(T=0\)。计算域长为\(100\)，初始条件为\(T=0\)。计算的参数设置为热扩散率\(\alpha=0.25\)，\(\delta_x=1\)，\(\delta_t=1\)，\(k=1\)。我们可以看到，随着时间的推移，温度会慢慢的变成线性。

2.3. 一维扩散方程的D1Q3模型

2.3.1. D1Q3模型整理：

三个离散速度

\[ e_0 = 0, \quad e_1 = c, \quad e_2 = -c \]

其中\(c=\delta_x/\delta_t\), \(\delta_x\)是空间步长，\(\delta_t\)是时间步长。

对应的权重

\[ w_0 = 4/6, \quad w_1 = w_2 = 1/6 \]

平衡态分布函数

\[ f_i^\mathrm{eq} = w_i T \]

演化方程

\[ f_i(\boldsymbol{x}+\boldsymbol{e}_i \delta_t,t+\delta_t)=f_i(\boldsymbol{x},t)-\frac{1}{\tau}[f_i(\boldsymbol{x},t)-f_i^{\mathrm{eq}}(\boldsymbol{x},t)] \]

也可以写成如下的形式：

\[ f_i(\boldsymbol{x}+\boldsymbol{e}_i \delta_t,t+\delta_t)=(1-\omega) f_i(\boldsymbol{x},t)+\omega f_i^{\mathrm{eq}}(\boldsymbol{x},t) \]

其中\(\omega = 1/\tau\)。

宏观量

\[ T= \sum_i f_i \]

热扩散率

\[ \alpha = \frac{\delta _x^2}{3\delta_t}\left(\tau - \frac{1}{2}\right) \]

2.3.2. 边界条件

边界条件是求解扩散方程的关键。在一维情况下，我们一般会有两种边界条件：

2.3.2.1. 第一类边界条件：恒温边界条件

这种边界条件是指在某个位置，温度是恒定的。比如在\(x=0\)处，\(T=1\)，在\(x=L\)处，\(T=0\)。这种边界条件的处理方法是：

左壁面

\[ f_1 = T_\mathrm{left} - f_2 - f_0 \]

右壁面

\[ f_2 = T_\mathrm{right} - f_1 - f_0 \]

2.3.2.2. 第二类边界条件：恒热流边界条件

这种边界条件是指在某个位置，热流密度是恒定的。比如在\(x=0\)处，$\(\frac{\partial T}{\partial x} = q\)$ 对于这种边界条件，我们需要知道的是，在D1Q2中，每个地方的热流密度可以通过以下方式计算：

\[ q = 3k\omega\frac{f_1 - f_2}{\delta_x} \]

其中，\(k\)为热导率。因此，我们可以通过以下方式来处理边界条件：

左壁面

\[ f_1 = f_2 + \frac{q\delta_x}{3k\omega} \]

右壁面

\[ f_2 = f_1 - \frac{q\delta_x}{3k\omega} \]

2.3.3. 算例验证

2.3.3.1. 恒温边界条件

左侧为\(T=1\)，右侧为\(T=0\)。计算域长为\(100\)，初始条件为\(T=0.5\)。计算的参数设置为热扩散率\(\alpha=0.25\)，\(\delta_x=1\)，\(\delta_t=1\)。我们可以看到，随着时间的推移，温度会慢慢的变成线性。

绘制热流密度的分布，FDM表示有限差分，计算方法为：

\[ q = -k \frac{\partial T}{\partial x} = -k \frac{T_{i+1}-T_{i}}{\delta_x} \]

LBM的计算方法为：

\[ q = 3k\omega \frac{f_1 - f_2}{\delta_x} \]

计算时，我们取\(k=1\)。可以看到两种方式计算的结果是一致的，而且与D1Q2的结果也是一致的。

2.3.3.2. 恒热流边界条件

此算例的边界条件为左侧恒定热流边界条件：\(q=0.01\)，右侧为恒定壁温边界条件：\(T=0\)。计算域长为\(100\)，初始条件为\(T=0\)。计算的参数设置为热扩散率\(\alpha=0.25\)，\(\delta_x=1\)，\(\delta_t=1\)，\(k=1\)。我们可以看到，随着时间的推移，温度会慢慢的变成线性。与D1Q2的结果也是一致的。

2.3.3.3. 变化的热扩散率

其实，前面的分析仅针对恒定热扩散率的情况，热扩散率时，上面介绍的LB模型所对应的宏观方程为：

\[ \frac{\partial T}{\partial t} = \frac{\partial}{\partial x}\left(\alpha \frac{\partial T}{\partial x}\right) \]

当\(\alpha\)为常数时，可以提出来，就变成了上述的形式。
使用LBM对变化的热扩散率求解时，只需要将\(\alpha\)替换成\(\alpha(x)\)即可。下面我们给出一个例子，\(\alpha(x)\)是一个线性函数：

\[ \alpha(x) = 0.1 + 0.009x \]

左侧为\(T=1\)，右侧为\(T=0\)。计算域长为\(100\)，初始条件为\(T=0.5\)。计算的参数设置为\(\delta_x=1\)，\(\delta_t=1\)。使用D1Q3模型求解，我们可以看到，随着时间的推移，左侧温度扩散的慢，右侧快。稳态时，斜率不一致，但热流密度\(q=\alpha \frac{\partial T}{\partial x}\)（有限差分方法计算得到）是相等的。 ​

2.3.3.4. 源项的处理

所谓源项，也就是在扩散方程中，右侧有一个非零的项。比如（热扩散率不变）：

\[ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} + S \]

这个方程使用LBM求解时，只需要在碰撞项中加上源项即可。源项的处理方法是：

\[ f_i(\boldsymbol{x}+\boldsymbol{e}_i \delta_t,t+\delta_t)=(1-\omega) f_i(\boldsymbol{x},t)+\omega f_i^{\mathrm{eq}}(\boldsymbol{x},t) + \delta_t\omega S \]

其他的地方都一致。下面我们给出一个例子，源项为：

\[ S = 0.0001 \]

两侧为\(T=0\)。计算域长为\(100\)，初始条件为\(T=0\)。计算的参数设置为热扩散率\(\alpha=0.25\)，\(\delta_x=1\)，\(\delta_t=1\)。使用D1Q3模型求解，我们可以看到，随着时间的推移，温度最后区域稳定。

使用有限差分计算\(q = \frac{\partial T}{\partial x}\)，认为\(k=1\)。

2.4. 模型validation和时空二阶精度检验

上面的计算只是一些定性的验证，我们有必要通过解析解来验证我们的模型。对于一维的扩散方程：

\[ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} \]

初始温度\(T=0\)，随后的时间左侧为\(T=0.0\)，右侧为\(T=1.0\)，\(L=100\),\(\alpha=0.25\)。我们可以得到解析解：

\[ T(x,t) = T_\mathrm{right}\frac{x}{L}-\frac{2T_\mathrm{right}}{\pi}\sum_{n=1}^{\infty}\frac{1}{n}\exp\left({-n^2\pi^2\frac{\alpha t}{L^2}}\right)\sin \left[n\pi\left(1-\frac{x}{L}\right)\right] \]

编程计算后，可以看到我们的计算结果与解析解十分吻合。

LBM所谓的具有时空二阶精度，是这么个情况，这里我们使用如下的公式来计算误差：

\[ \epsilon = \sqrt{\frac{1}{N}\sum_{i=1}^{N}\left(T_{\mathrm{ana}}(x_i,t) - T_{\mathrm{num}}(x_i,t)\right)^2} \]

其中\(T_{\mathrm{ana}}\)是解析解，\(T_{\mathrm{num}}\)是数值解。我们可以看到，随着\(\delta_x\)和\(\delta_t\)的减小，误差是会减小的。
我们固定\(\omega=0.8\)，则\(\delta_t\)与\(\delta_x\)的关系为：

\[ \delta_t = \frac{\delta_x^2}{3\alpha}\left(\tau - \frac{1}{2}\right)=\frac{\delta_x^2}{3\alpha}\left(\frac{1}{\omega} - \frac{1}{2}\right) \]

选取不同的\(\delta_x\)和\(\delta_t,t=5000\)时的误差)：

\(\delta_x\)	\(\delta_t\)	\(\epsilon\)
10	100	1.81e-03
1	1	1.95e-05
0.1	0.01	1.96e-07

我们可以看到，随着\(\delta_x\)和\(\delta_t\)的减小，误差是会减小的。

画出\(\log(\epsilon)\)与\(\log(\delta_x)\)的关系，我们可以看到，误差与\(\delta_x\)的关系是线性的，斜率十分接近2，这就表明我们的模型是时空二阶精度的。

2.5. 二维扩散模型的D2Q5模型

二维的扩散方程是这样的：

\[ \frac{\partial T}{\partial t} = \alpha \nabla^2 T \]

写成离散形式：

\[ \frac{\partial T}{\partial t} = \alpha \left(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}\right) \]

2.5.1. 模型整理

我们可以使用D2Q5模型来求解这个方程。D2Q5模型的整理如下：

五个离散速度

\[\begin{split} e_i = \left\{ \begin{aligned} &0, \quad i=0 \\ &c\left[\cos \frac{\pi (i-1)}{2}, \sin \frac{\pi (i-1)}{2}\right], \quad i=1,2,3,4 \end{aligned} \right. \end{split}\]

对应的权重

\[ w_0 = 2/6, \quad w_{1-4} = 1/6 \]

平衡态分布函数

\[ f_i^\mathrm{eq} = w_i T \]

演化方程

\[ f_i(\boldsymbol{x}+\boldsymbol{e}_i \delta_t,t+\delta_t)=f_i(\boldsymbol{x},t)-\frac{1}{\tau}[f_i(\boldsymbol{x},t)-f_i^{\mathrm{eq}}(\boldsymbol{x},t)] \]

也可以写成如下的形式：

\[ f_i(\boldsymbol{x}+\boldsymbol{e}_i \delta_t,t+\delta_t)=(1-\omega) f_i(\boldsymbol{x},t)+\omega f_i^{\mathrm{eq}}(\boldsymbol{x},t) \]

其中\(\omega = 1/\tau\)。

宏观量

\[ T= \sum_i f_i \]

热扩散率

\[ \alpha = \frac{\delta _x^2}{3\delta_t}\left(\tau - \frac{1}{2}\right) \]

2.5.2. 边界条件

边界条件是求解扩散方程的关键。这里主要介绍恒温边界条件，恒热流边界条件可以转化为恒温边界条件。 壁温为\(T_\mathrm{w}\)，我们可以通过以下方式处理：

左壁面（\(f_1\)未知）

\[ f_1 = T_\mathrm{left} - \sum_{i=0,2,3,4}f_i \]

右壁面（\(f_3\)未知）

\[ f_3 = T_\mathrm{right} - \sum_{i=0,1,2,4}f_i \]

上壁面（\(f_4\)未知）

\[ f_4 = T_\mathrm{top} - \sum_{i=0,1,2,3}f_i \]

下壁面（\(f_2\)未知）

\[ f_2 = T_\mathrm{bottom} - \sum_{i=0,1,3,4}f_i \]

2.5.3. 算例验证

本算例为一个\(50\times50\)的计算域，初始时刻温度为0，左壁面温度为1，右壁面温度为0，上壁面温度为0，下壁面为绝热。算例的参数设置为\(\delta_x=1\)，\(\delta_t=1\)，\(\alpha=0.1​\)，计算至\(t=2000​\)。

该模型\(\delta_t\)和\(\delta_x\)均可以任意调节。

2.6. 二维扩散方程的D2Q9模型

2.6.1. 模型整理

我们可以使用D2Q9模型来求解这个方程。D2Q9模型的整理如下：

九个离散速度

\[\begin{split} e_i = \left\{ \begin{aligned} &0, \quad i=0 \\ &c\left[\cos \frac{\pi (i-1)}{2}, \sin \frac{\pi (i-1)}{2}\right], \quad i=1,2,3,4 \\ &\sqrt{2}c\left( \cos \left[ (2i -1)\frac{\pi}{4} \right] ,\sin \left[ (2i -1)\frac{\pi}{4} \right] \right), \quad i=5,6,7,8 \end{aligned} \right. \end{split}\]

对应的权重

\[ w_0 = 4/9, \quad w_{1-4} = 1/9, \quad w_{5-8} = 1/36 \]

平衡态分布函数

\[ f_i^\mathrm{eq} = w_i T \]

演化方程

\[ f_i(\boldsymbol{x}+\boldsymbol{e}_i \delta_t,t+\delta_t)=f_i(\boldsymbol{x},t)-\frac{1}{\tau}[f_i(\boldsymbol{x},t)-f_i^{\mathrm{eq}}(\boldsymbol{x},t)] \]

也可以写成如下的形式：

\[ f_i(\boldsymbol{x}+\boldsymbol{e}_i \delta_t,t+\delta_t)=(1-\omega) f_i(\boldsymbol{x},t)+\omega f_i^{\mathrm{eq}}(\boldsymbol{x},t) \]

其中\(\omega = 1/\tau\)。

宏观量

\[ T= \sum_i f_i \]

热扩散率

\[ \alpha = \frac{\delta _x^2}{3\delta_t}\left(\tau - \frac{1}{2}\right) \]

2.6.2. 边界条件

这里的边界条件就比较困难了，因为对于一个壁面，未知量有三个，例如左壁面，\(f_1\)，\(f_5\)，\(f_8\)都是未知的。但我们只知道壁面的温度。对于温度场，非平衡态外推格式是一个非常好用的格式。该格式的基本假设为，壁面处的非平衡态部分和临近格子的非平衡态部分是一致的，即：

\[ f_\mathrm{near} - f^\mathrm{eq}_\mathrm{near} = f_\mathrm{wall} - f^\mathrm{eq}_\mathrm{wall} \]

其中，near表示临界格子，wall表示壁面格子。则：

\[ f_\mathrm{wall} = f^\mathrm{eq}_\mathrm{wall} + f_\mathrm{near} - f^\mathrm{eq}_\mathrm{near} \]

2.6.3. 算例验证

本算例与5.3节一致，只是将D2Q5模型改为D2Q9模型。计算结果如下：

2.7. 总结

本节我们讲解了扩散方程的求解，包括一维和二维的情况。我们使用了D1Q2，D1Q3，D2Q5，D2Q9模型来求解扩散方程。并且通过解析解验证了我们的模型。最后我们还验证了我们的模型是时空二阶精度的。希望大家通过这节课的学习，对LBM的扩散方程有了一个直观的认识。

2.8. 代码下载

扩散方程代码.zip

---

如果您觉得这个项目对您有帮助，可以考虑用以下方式支持我：

• ☕ 请我喝杯咖啡

如果条件允许，欢迎捐赠支持！
每一分都是对我莫大的鼓励，让我能投入更多时间维护和更新。

• 📄 引用我的文章

如果暂时不便捐赠，适当引用我的文章也是极好的支持！
您的引用能帮助这个工作获得更多关注，同样让我感到无比欣慰。

无论哪种方式，都是对我的巨大支持！🙏感谢您让开源世界更美好！✨