# 扩散方程 文档B站视频讲解： > [LBM入门课程(2.1)LBM求解扩散方程(一维)](https://www.bilibili.com/video/BV1oc2HYyEgo/?share_source=copy_web&vd_source=6cb794fb9d90103946fb79ba5f60451f) > [LBM入门课程(2.2)LBM的时空二阶精度及变系数/带源项的扩散方程求解](https://www.bilibili.com/video/BV1EimNY6E7w/?share_source=copy_web&vd_source=6cb794fb9d90103946fb79ba5f60451f) > [LBM入门课程(2.3)D2Q5/D2Q9求扩散方程及非平衡态外推边界条件](https://www.bilibili.com/video/BV1pvSHYAER5/?share_source=copy_web&vd_source=6cb794fb9d90103946fb79ba5f60451f) 这节我们讲解**扩散方程**的格子玻尔兹曼方法的实现。包括一维和二维的情况。 ## 什么是扩散方程？ 扩散方程是描述物质扩散的方程。在一维情况下，扩散方程是这样的： $$ \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$的边界条件，然后就会慢慢的向右边扩散。这就是称它为扩散方程的原因。 ```{figure} fig/fig2/扩散方程图片1.png :width: 50% :alt: logo :align: center ``` 下图则展示了一个初始条件为正弦分布场随时间的演化，两侧为周期性边界条件。 ```{figure} fig/fig2/扩散方程图片2.png :width: 50% :alt: logo :align: center ``` 相信大家通过这两个算例已经对扩散方程有了一个直观的认识。接下来我们就来看看如何用LBM来求解扩散方程的。 ## 一维扩散方程的D1Q2模型 重新写一个一维扩散方程： $$ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} $$ 当把$T$当作变量时，$\alpha$一般称为热扩散系数。 ### 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) $$ ### 边界条件 边界条件是求解扩散方程的关键。在一维情况下，我们一般会有两种边界条件： #### 第一类边界条件：恒温边界条件 这种边界条件是指在某个位置，温度是恒定的。比如在$x=0$处，$T=1$，在$x=L$处，$T=0$。这种边界条件的处理方法是： > 左壁面 $$ f_1 = T_\mathrm{left} - f_2 $$ > 右壁面 $$ f_2 = T_\mathrm{right} - f_1 $$ #### 第二类边界条件：恒热流边界条件 这种边界条件是指在某个位置，热流密度是恒定的。比如在$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} $$ ### 算例验证 #### 恒温边界条件 左侧为$T=1$，右侧为$T=0$。计算域长为$100$，初始条件为$T=0.5$。计算的参数设置为热扩散率$\alpha=0.25$，$\delta_x=1$，$\delta_t=1$。我们可以看到，随着时间的推移，温度会慢慢的变成线性。 ```{figure} fig/fig2/D1Q2_T_2.png :width: 50% :alt: logo :align: center ``` 绘制热流密度的分布，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$。可以看到两种方式计算的结果是一致的。 ```{figure} fig/fig2/D1Q2_T_1.png :width: 50% :alt: logo :align: center ``` 下面展示的是$\delta_x = 0.5$，$\delta_t=0.5$时的计算结果： ```{figure} fig/fig2/D1Q2_T_3.png :width: 50% :alt: logo :align: center ``` ```{figure} fig/fig2/D1Q2_T_4.png :width: 50% :alt: logo :align: center ``` 可以看到我们的程序可以独立的调节$\delta_x$和$\delta_t$，并且计算结果是正确的。 #### 恒热流边界条件 此算例的边界条件为左侧恒定热流边界条件：$q=0.01$，右侧为恒定壁温边界条件：$T=0$。计算域长为$100$，初始条件为$T=0$。计算的参数设置为热扩散率$\alpha=0.25$，$\delta_x=1$，$\delta_t=1$，$k=1$。我们可以看到，随着时间的推移，温度会慢慢的变成线性。 ```{figure} fig/fig2/D1Q2_q_1.png :width: 50% :alt: logo :align: center ``` ```{figure} fig/fig2/D1Q2_q_2.png :width: 50% :alt: logo :align: center ``` ## 一维扩散方程的D1Q3模型 ### 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) $$ ### 边界条件 边界条件是求解扩散方程的关键。在一维情况下，我们一般会有两种边界条件： #### 第一类边界条件：恒温边界条件 这种边界条件是指在某个位置，温度是恒定的。比如在$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 $$ #### 第二类边界条件：恒热流边界条件 这种边界条件是指在某个位置，热流密度是恒定的。比如在$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} $$ ### 算例验证 #### 恒温边界条件 左侧为$T=1$，右侧为$T=0$。计算域长为$100$，初始条件为$T=0.5$。计算的参数设置为热扩散率$\alpha=0.25$，$\delta_x=1$，$\delta_t=1$。我们可以看到，随着时间的推移，温度会慢慢的变成线性。 ```{figure} fig/fig2/D1Q3_T_1.png :width: 50% :alt: logo :align: center ``` 绘制热流密度的分布，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的结果也是一致的。 ```{figure} fig/fig2/D1Q3_T_2.png :width: 50% :alt: logo :align: center ``` #### 恒热流边界条件 此算例的边界条件为左侧恒定热流边界条件：$q=0.01$，右侧为恒定壁温边界条件：$T=0$。计算域长为$100$，初始条件为$T=0$。计算的参数设置为热扩散率$\alpha=0.25$，$\delta_x=1$，$\delta_t=1$，$k=1$。我们可以看到，随着时间的推移，温度会慢慢的变成线性。与D1Q2的结果也是一致的。 ```{figure} fig/fig2/D1Q3_q_1.png :width: 50% :alt: logo :align: center ``` ```{figure} fig/fig2/D1Q3_q_2.png :width: 50% :alt: logo :align: center ``` #### 变化的热扩散率 其实，前面的分析仅针对恒定热扩散率的情况，热扩散率时，上面介绍的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}$（有限差分方法计算得到）是相等的。 ​ ```{figure} fig/fig2/D1Q3_alpha_1.png :width: 50% :alt: logo :align: center ``` ```{figure} fig/fig2/D1Q3_alpha_2.png :width: 50% :alt: logo :align: center ``` #### 源项的处理 所谓源项，也就是在扩散方程中，右侧有一个非零的项。比如（热扩散率不变）： $$ \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模型求解，我们可以看到，随着时间的推移，温度最后区域稳定。 ```{figure} fig/fig2/D1Q3_S_1.png :width: 50% :alt: logo :align: center ``` ```{figure} fig/fig2/D1Q3_S_2.png :width: 50% :alt: logo :align: center ``` 使用有限差分计算$q = \frac{\partial T}{\partial x}$，认为$k=1$。 ## 模型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] $$ 编程计算后，可以看到我们的计算结果与解析解十分吻合。 ```{figure} fig/fig2/D1Q3_validation.png :width: 50% :alt: logo :align: center ``` 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$的减小，误差是会减小的。 ```{figure} fig/fig2/D1Q3_error_dx_10.png :width: 50% :alt: logo :align: center ``` ```{figure} fig/fig2/D1Q3_error_dx_1.png :width: 50% :alt: logo :align: center ``` ```{figure} fig/fig2/D1Q3_error_dx_0.1.png :width: 50% :alt: logo :align: center ``` 画出$\log(\epsilon)$与$\log(\delta_x)$的关系，我们可以看到，误差与$\delta_x$的关系是线性的，斜率十分接近2，这就表明我们的模型是**时空二阶精度**的。 ```{figure} fig/fig2/D1Q3_log-log.png :width: 50% :alt: logo :align: center ``` ## 二维扩散模型的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) $$ ### 模型整理 我们可以使用D2Q5模型来求解这个方程。D2Q5模型的整理如下： > 五个离散速度 $$ 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. $$ > 对应的权重 $$ 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) $$ ### 边界条件 边界条件是求解扩散方程的关键。这里主要介绍恒温边界条件，恒热流边界条件可以转化为恒温边界条件。 壁温为$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 $$ ### 算例验证 本算例为一个$50\times50$的计算域，初始时刻温度为0，左壁面温度为1，右壁面温度为0，上壁面温度为0，下壁面为绝热。算例的参数设置为$\delta_x=1$，$\delta_t=1$，$\alpha=0.1​$，计算至$t=2000​$。 ```{figure} fig/fig2/D2Q5.png :width: 50% :alt: logo :align: center ``` 该模型$\delta_t$和$\delta_x$均可以任意调节。 ## 二维扩散方程的D2Q9模型 ### 模型整理 我们可以使用D2Q9模型来求解这个方程。D2Q9模型的整理如下： > 九个离散速度 $$ 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. $$ > 对应的权重 $$ 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) $$ ### 边界条件 这里的边界条件就比较困难了，因为对于一个壁面，未知量有三个，例如左壁面，$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} $$ ### 算例验证 本算例与5.3节一致，只是将D2Q5模型改为D2Q9模型。计算结果如下： ```{figure} fig/fig2/D2Q9.png :width: 50% :alt: logo :align: center ``` ## 总结 本节我们讲解了扩散方程的求解，包括一维和二维的情况。我们使用了D1Q2，D1Q3，D2Q5，D2Q9模型来求解扩散方程。并且通过解析解验证了我们的模型。最后我们还验证了我们的模型是时空二阶精度的。希望大家通过这节课的学习，对LBM的扩散方程有了一个直观的认识。 ## 代码下载 [扩散方程代码.zip](https://raw.githubusercontent.com/JR-LBM-THU/JRLBM-tutorial/main/docs/source/LBM基础入门/files/2.扩散方程代码.zip) --- **如果您觉得这个项目对您有帮助，可以考虑用以下方式支持我：** - ☕ 请我喝杯咖啡 > 如果条件允许，[欢迎捐赠支持](../支持与捐赠/支持与捐赠.md)！ > 每一分都是对我莫大的鼓励，让我能投入更多时间维护和更新。 - 📄 引用我的文章 > 如果暂时不便捐赠，[适当引用我的文章](../我的文章/我的文章.md)也是极好的支持！ > 您的引用能帮助这个工作获得更多关注，同样让我感到无比欣慰。 无论哪种方式，都是对我的巨大支持！🙏感谢您让开源世界更美好！✨