3. 对流扩散方程

文档B站视频讲解：

格子玻尔兹曼方法(LBM)求解对流扩散方程

这节我们介绍一下对流扩散方程的格子玻尔兹曼方法实现。包括一维和二维。

3.1. 什么是对流扩散方程？

扩散方程我们之前已经进行了介绍。那我们首先来看一下什么是对流方程。对于一维问题，对流方程是这样的：

\[ \frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x}=0 \]

其中\(u\)表示对流速度。我们知道\(\phi(x,t)\)是空间和时间的函数，所以上式其实表征的意思是\(\phi\)对时间的全导数，即:

\[ \frac{\mathrm{d} \phi(x,t)}{\mathrm{d} t}=\frac{\partial \phi}{\partial t} + \frac{\partial x}{\partial t} \frac{\partial \phi}{\partial x}=\frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x}=0 \]

这里\(u=\partial x/\partial t\)，表征\(\phi\)移动的速度。下图是一个对流方程结果的经典展示：

图中\(u=1\)

那么对流扩散方程，就是既包含对流项，也包含扩散。一维情况下，方程如下：

\[ \frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} = \alpha \frac{\partial^2 \phi}{\partial x^2} \]

写成高维矢量形式维：

\[ \frac{\mathrm{d} \phi(\boldsymbol{x},t)}{\mathrm{d} t}=\frac{\partial \phi}{\partial t} +\boldsymbol{u}\cdot \nabla \phi = \alpha \nabla^2 \phi \]

这里的\(\boldsymbol{x}\)表示位置矢量，\(\boldsymbol{u}\)表示速度矢量。展开到二维：

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

注，其实将\(\phi\)换成\(T\)就是质量、动量、能量方程中的能量方程。将\(\phi\)换成\(\rho \boldsymbol{u}\)就是动量方程。

3.2. 一维对流扩散方程的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 \]

平衡态分布函数（\(c_s^2=c^2/3\)）

\[ f_i^\mathrm{eq} = w_i \phi\left( 1 + \frac{e_i\cdot u}{c^2_s}\right) \]

演化方程

\[ 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\)。

宏观量

\[ \phi= \sum_i f_i \]

扩散系数

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

3.2.1. 算例1模拟结果

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

图中\(u=0.1\)

图中\(u=0.05\)

图中\(u=0.0\)

3.2.2. 算例2模拟结果

两侧周期性边界条件。计算域长为200，\(10<x<30\)内\(\phi=1\)，其他区域\(\phi=0\)。计算的参数设置为热扩散率\(\alpha=0.05\)，\(\delta_x=1\)，\(\delta_t=1\)。我们可以看到，随着时间的推移，温度会慢慢的变化。

图中\(u=0.1\)

图中\(u=0.05\)

图中\(u=0.0\)

3.3. 二维对流扩散方程的D2Q5模型

从上面的与扩散方程的比较看出，使用D1Q3模型求解对流扩散方程与求解扩散方程唯一的不同就是平衡态分布函数。同样的，D2Q5也如此，平衡态分布函数为：

\[ f_i^\mathrm{eq} = w_i \phi\left( 1 + \frac{\boldsymbol{e_i\cdot u}}{c^2_s}\right) \]

其他的内容与上一章介绍的D2Q5求解扩散方程一致。

3.3.1. 模拟结果

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

图中\(u_x=0.0\)

图中\(u_x=0.02\)

图中\(u_x=0.04\)

图中\(u_x=0.06\)

图中\(u_x=0.08\)

图中\(u_x=0.10\)

比较不同时刻的结果，对流与扩散的效果清晰可见~

3.4. 二维对流扩散方程的D2Q9模型

从上面的与扩散方程的比较看出，使用D1Q3模型求解对流扩散方程与求解扩散方程唯一的不同就是平衡态分布函数。同样的，D2Q9也如此，平衡态分布函数为：

\[ f_i^\mathrm{eq} = w_i \phi\left( 1 + \frac{\boldsymbol{e_i\cdot u}}{c^2_s}\right) \]

其他的内容与上一章介绍的D2Q9求解扩散方程一致。

3.4.1. 模拟结果

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

图中\(u_x=0.0\)

图中\(u_x=0.02\)

图中\(u_x=0.04\)

图中\(u_x=0.06\)

图中\(u_x=0.08\)

图中\(u_x=0.10\)

比较不同时刻的结果，对流与扩散的效果清晰可见~比较D2Q9和D2Q5的结果，可以看到两个模型计算的结果也都是一致的。

\(u_x=0.05\)时下边界的温度值

3.5. 代码下载

对流扩散代码.zip

---

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

• ☕ 请我喝杯咖啡

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

• 📄 引用我的文章

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

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