💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

🎁完整资源、论文复现、期刊合作、论文辅导及科研仿真定制事宜点击：

👉👉👉本文完整资源下载

⛳️座右铭：行百里者，半于九十。

⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能解答你胸中升起的一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

💥第一部分——内容介绍

基于格子Boltzmann方法的二维通道多孔介质流动模拟研究

摘要

多孔介质内的流体输运行为在油气开采、地下水治理、化学反应工程以及生物组织灌注等诸多领域具有重要的科学价值与工程意义。受限于复杂孔隙几何带来的网格生成难题，传统连续介质数值方法在孔隙尺度模拟中往往面临效率与精度的双重挑战。本文以介观尺度的格子Boltzmann方法为核心工具，构建了一套面向二维通道内多孔介质流动的数值模拟框架。研究采用九速度离散模型描述流体粒子的迁移与碰撞过程，借助Bhatnagar–Gross–Krook近似简化碰撞算子，并分别引入反弹格式处理无滑移固体壁面与浸入障碍物，引入Zou–He格式施加进出口的速度与压力约束。在此基础上，本文依次模拟了平板通道内的Poiseuille流动、圆柱与方柱绕流、真实裂缝图像驱动的多孔介质渗流以及正弦孔道流动，并将模拟结果与解析解和公开文献基准数据进行系统对比。收敛性分析表明，所采用的数值格式在空间上具有接近二阶的精度；稳定性分析则揭示了马赫数、松弛时间与雷诺数之间的相互制约关系。研究结果表明，格子Boltzmann方法与浸入边界思想的结合，为复杂几何条件下的孔隙尺度流动提供了一条算法简洁、边界灵活、易于并行的高效求解路径。

关键词：格子Boltzmann方法；多孔介质；浸入边界；反弹边界；Zou–He边界；孔隙尺度流动

1. 引言

1.1 研究背景

多孔介质是由固体骨架与相互连通的孔隙空间共同构成的多相体系，其内部流体流动广泛存在于自然界与工业过程中。在石油工程领域，原油、天然气与地层水在岩石孔隙中的渗流规律直接决定了油藏的采收率与开发方案；在水文地质领域，污染物在含水层中的迁移扩散依赖于孔隙结构的连通性；在化学工程领域，多孔催化剂与填充床反应器的传质效率受孔隙尺度流场分布的显著影响；在生物力学领域，血液在毛细血管网络与骨组织中的流动同样可视为多孔介质渗流问题。因此，准确刻画孔隙尺度下的速度场与压力场，对于理解宏观输运规律、建立可靠的升尺度模型具有基础性意义。

1.2 传统方法的局限

传统计算流体力学方法在模拟多孔介质流动时主要面临三方面困难。其一，孔隙几何高度不规则，直接生成贴体网格需要耗费大量人力与计算资源；其二，孔隙尺度与宏观尺度跨越多个数量级，全域直接模拟的计算代价极为高昂；其三，固体壁面与流体域的交界处往往存在复杂的局部流动结构，对数值格式的精度与稳定性提出了苛刻要求。尽管有限体积法与有限元法在规则几何中表现优异，但在处理随机孔隙结构时仍显笨拙。

1.3 格子Boltzmann方法的优势

格子Boltzmann方法作为一种介于微观分子动力学与宏观连续介质力学之间的介观数值方法，自二十世纪九十年代以来迅速发展。该方法以分布函数的演化代替对宏观守恒方程的直接离散，具有算法结构简单、天然适合并行计算、边界处理灵活等突出优点。尤其值得指出的是，该方法通过简单的反弹规则即可实现无滑移边界，通过布尔掩膜即可描述任意复杂的固体几何，这为多孔介质流动模拟提供了极大的便利。近年来，已有大量研究将格子Boltzmann方法应用于岩石孔隙渗流、裂缝流动、相对渗透率预测以及多相驱替等问题，取得了丰硕成果。

1.4 本文工作

本文在前人研究基础上，系统构建了一套基于格子Boltzmann方法的二维通道多孔介质流动模拟框架。研究按照由简到繁的逻辑逐步推进：首先验证平板通道内的Poiseuille流动，确保基础求解器正确；随后引入浸入边界条件，模拟圆柱与方柱绕流；进而导入真实裂缝图像，模拟多孔介质渗流；最后构造正弦孔道，考察变截面通道内的流动特征。全文重点关注数值格式的精度、稳定性与适用边界，并对模拟结果进行多角度验证。

2. 格子Boltzmann方法理论基础

2.1 从分子运动到介观分布函数

流体在微观尺度上由大量离散分子组成，分子的随机运动与碰撞决定了宏观流动行为。直接追踪每一个分子的运动轨迹既不现实也无必要。Boltzmann方程通过对分子分布函数进行统计描述，在保留分子运动统计特征的同时大幅降低了自由度。分布函数表示在某一时刻、某一位置附近、某一速度区间内找到分子的概率密度。通过对分布函数求矩，可以得到密度、动量与能量等宏观物理量。

格子Boltzmann方法在Boltzmann方程的基础上进一步引入两重离散：速度空间被离散为有限个方向，物理空间被离散为规则格点，时间被离散为等间隔步长。经过这样的离散化，原本连续的分布函数演化方程转化为一组代数递推关系，便于计算机实现。

2.2 碰撞与迁移的双步机制

格子Boltzmann方法的核心演化过程可以形象地理解为“碰撞—迁移”两个交替步骤。在碰撞步骤中，位于同一格点上的不同方向分布函数相互交换动量，使系统向局部平衡态趋近；在迁移步骤中，碰撞后的分布函数沿各自对应的离散速度方向移动到相邻格点。这一机制既体现了分子碰撞的物理本质，又保证了质量与动量的局部守恒。

碰撞过程的建模方式决定了方法的精度与稳定性。本文采用最经典的BGK近似，即假设碰撞后分布函数以单一松弛时间线性趋近于平衡态。该近似形式简洁、物理意义明确，适用于低马赫数与低雷诺数范围内的不可压缩流动。平衡态分布函数本身由局部密度与速度决定，其具体形式通过对麦克斯韦–玻尔兹曼分布进行低马赫数泰勒展开得到。

2.3 离散速度模型的选择

在二维情形下，最常用的离散速度模型为九速度模型，即每个格点包含一个静止方向、四个轴向方向与四个对角方向。这一模型在保证足够对称性的同时，能够正确恢复二维不可压缩Navier–Stokes方程。每个方向对应一个权重系数，权重的大小反映了该方向对平衡态分布的贡献程度。静止方向的权重最大，轴向方向次之，对角方向最小。

2.4 宏观量的提取

在每个时间步，通过对所有方向的分布函数求和即可获得局部密度；通过将各方向的分布函数与其速度向量相乘后求和，再除以密度，即可获得局部速度。这一过程本质上是分布函数的零阶与一阶矩运算，物理意义清晰，计算实现简便。

2.5 输运系数与松弛时间

通过多尺度展开分析可以证明，格子Boltzmann方法在低马赫数极限下能够恢复不可压缩Navier–Stokes方程，且运动黏度与松弛时间之间存在确定的代数关系。当松弛时间趋近于二分之一时，黏度趋近于零，对应高雷诺数流动；当松弛时间远大于二分之一时，黏度较大，对应低雷诺数流动。为保证数值稳定，松弛时间必须严格大于二分之一。这一约束在高雷诺数模拟中尤为关键，因为此时松弛时间接近临界值，碰撞过程变慢，数值振荡容易累积并导致发散。

3. 边界条件处理

3.1 反弹边界与无滑移条件

在固体壁面处，流体速度必须满足无滑移条件，即壁面处流体速度等于壁面速度。在格子Boltzmann框架下，实现无滑移最简单的方式是反弹规则：当分布函数迁移至边界格点时，沿原方向的反方向反射回流体域。对于静止壁面，反射后的分布函数等于入射方向的反方向分布函数。这一规则物理上对应于流体粒子与壁面碰撞后原路返回，宏观上自然满足无滑移条件。

反弹边界可分为格点反弹与半格点反弹两种实现方式。前者将边界置于格点之上，实现简单但仅具有一阶精度；后者将边界置于格点之间，具有二阶精度但实现稍复杂。本文在多孔介质模拟中采用格点反弹，主要考虑其便于与图像掩膜结合的优点。

3.2 浸入边界与布尔掩膜

浸入边界方法的基本思想是将复杂固体几何嵌入规则网格之中，通过标记固体格点与流体格点来区分不同区域。在格子Boltzmann框架下，只需构造一个与网格同尺寸的布尔矩阵，将固体区域标记为一，流体区域标记为零，即可实现任意复杂几何的描述。在碰撞步骤之后，对所有固体格点施加反弹规则，即可实现无滑移约束。这一方法的突出优势在于：无需生成贴体网格，几何描述与数值求解完全解耦，特别适合处理随机孔隙、裂缝网络与生物组织等复杂结构。

3.3 Zou–He速度与压力边界

在通道流动模拟中，入口与出口边界条件的选取直接影响质量守恒与数值稳定性。本文采用Zou–He边界条件，在入口施加速度剖面，在出口施加压力约束。入口速度剖面取为抛物线分布，满足通道中心速度最大、壁面速度为零的Poiseuille流特征。出口给定恒定密度，相当于给定压力差，驱动流体持续流动。

Zou–He边界条件的核心思想是利用已知的宏观量，结合分布函数的矩关系，反推未知方向的分布函数。对于入口边界，已知水平速度与垂直速度，未知密度与部分方向分布函数，通过求解线性方程组即可闭合。对于出口边界，已知密度，未知速度，同样可以通过矩关系求解。采用压力边界而非速度边界作为出口条件，可以有效避免质量不守恒导致的速度漂移，有利于长时间模拟的稳定性。

3.4 周期性边界与体积力驱动

作为Zou–He边界的替代方案，周期性边界条件将出口的分布函数直接复制到入口，适用于充分发展的流动。在周期性边界下，流动通常由体积力驱动。体积力在每一时间步均匀添加到所有流体格点，相当于在动量方程中引入恒定压力梯度。体积力的大小由目标平均速度、运动黏度与通道宽度共同决定。本文以Zou–He边界为主，周期性边界作为备选方案，未做深入展开。

4. 几何建模与多孔介质构造

4.1 规则障碍物

为验证浸入边界方法的正确性，本文首先模拟了圆柱与方柱绕流。圆柱通过圆心坐标与半径定义，方柱通过中心坐标与边长定义。在布尔掩膜中，满足几何条件的格点被标记为固体，其余为流体。圆柱绕流是经典的流体力学基准问题，其尾涡结构与阻力系数已有大量文献数据可供对比。方柱绕流则因尖角处流动分离更为剧烈，对数值格式的稳定性提出了更高要求。

4.2 正弦孔道

为模拟岩石中常见的孔喉结构，本文构造了上下壁面均呈正弦变化的通道。正弦孔道的振幅与波数可调，能够模拟不同孔喉比与不同孔隙数量的情形。在孔喉处，流通面积减小，速度增大；在孔腹处，流通面积增大，速度减小。这一几何特征与真实岩石孔隙中的流动特征高度相似，可用于考察变截面通道内的流线分布与压力变化。

4.3 图像驱动的真实多孔介质

本文最具特色之处在于利用真实裂缝图像驱动多孔介质流动模拟。具体做法是：读取裂缝图像，将其转换为二值图像，其中裂缝区域为流体，岩石骨架为固体；然后将二值图像直接作为布尔掩膜输入求解器。这一流程无需任何几何重构，即可实现真实岩石裂缝网络的流动模拟。该方法的优势在于：可以直接利用实验获取的薄片图像、CT扫描切片或扫描电镜图像，实现从真实几何到数值模拟的无缝衔接。其局限性在于：图像分辨率决定了最小可分辨孔隙尺度，且二维切片无法完全反映三维连通性。

5. 数值实现流程

5.1 初始化

模拟开始时，所有格点的密度设为均匀值，速度场设为零或给定的初始剖面。分布函数初始化为对应局部密度与速度的平衡态分布。入口速度剖面按抛物线分布赋给入口边界格点。

5.2 迭代主循环

每一次迭代包含以下步骤：首先由分布函数求和计算密度与速度；然后施加入口与出口边界条件，更新边界格点的分布函数；接着计算平衡态分布函数，并执行碰撞操作，使分布函数向平衡态松弛；随后对浸入边界格点施加反弹规则；最后执行迁移操作，将分布函数沿各离散速度方向移动到相邻格点。上述过程不断重复，直至速度场变化小于预设阈值。

5.3 收敛判据

采用相邻两次迭代速度场的范数差作为收敛指标。当水平速度与垂直速度的变化量均小于给定容差时，认为流动达到稳态。容差的选取需在精度与计算成本之间权衡：容差过大会导致过早终止，结果未充分发展；容差过小则计算时间显著增加。

5.4 阻力系数计算

对于绕流问题，阻力系数是重要的验证指标。本文采用动量交换法计算流体对固体的作用力：对固体格点周围所有方向的分布函数与平衡态分布函数求和，并乘以对应方向的速度分量，即可得到水平与垂直方向的合力。阻力系数由合力除以动压与特征面积得到。计算过程中每隔若干时间步记录一次，取稳态后的峰值或平均值作为最终结果。

6. 结果分析与验证

6.1 Poiseuille流动验证

在无任何障碍物的平板通道中，Zou–He边界条件驱动流体形成抛物线速度剖面。将模拟得到的通道中心线上的水平速度与解析解对比，两者高度吻合。最大偏差出现在壁面附近，主要源于格点反弹边界的一阶精度。整体而言，基础求解器的正确性得到确认。

6.2 圆柱绕流

在中等雷诺数下，圆柱后方出现对称的分离泡；随着雷诺数增大，分离泡拉长并最终失稳，形成交替脱落的涡街。本文模拟得到的流线图清晰展示了尾涡结构，阻力系数随雷诺数的变化趋势与文献中的标准阻力曲线一致。在雷诺数约为一百五十时，模拟得到的阻力系数与文献值接近，验证了浸入边界方法的有效性。

6.3 方柱绕流

方柱绕流因尖角处流动分离更为剧烈，其尾涡结构比圆柱更为复杂。模拟结果显示，方柱后方出现周期性涡脱落，阻力系数随时间呈振荡变化。振荡频率对应的斯特劳哈尔数与文献值吻合。阻力系数的平均值与文献中的数值结果接近，进一步验证了方法在处理尖锐几何时的稳健性。

6.4 裂缝介质流动

导入真实裂缝图像后，模拟得到的速度场呈现出强烈的非均匀性。裂缝开度较大的区域速度较高，开度较小的区域速度较低；裂缝交汇处出现复杂的局部流动结构。与文献中基于同类图像的速度矢量图对比，整体分布趋势一致。这一结果说明，格子Boltzmann方法结合图像掩膜可以有效模拟真实岩石裂缝中的流动。

6.5 正弦孔道流动

在正弦孔道中，流线在孔喉处密集，在孔腹处稀疏，表明速度在孔喉处增大、在孔腹处减小。模拟得到的流线形态与文献结果整体一致，差异主要源于流线绘制方式的不同：文献中的流线从孔喉处起始，而本文的流线绘制方法仅追踪水平速度方向，导致部分流线未能完整展现回流结构。尽管如此，速度场的整体分布仍验证了方法的正确性。

7. 精度与稳定性分析

7.1 空间收敛性

为评估数值格式的空间精度，本文在三套不同分辨率的网格上分别求解同一问题，并计算相邻网格间误差的比值。结果表明，收敛阶接近二阶，与理论预期一致。若采用半格点反弹边界替代格点反弹，收敛阶有望进一步提高。格点反弹的一阶精度主要源于边界位置与格点位置的偏差，在半格点反弹中这一偏差被消除。

7.2 稳定性约束

格子Boltzmann方法的稳定性受多重因素制约。首先，马赫数必须远小于一，以抑制可压缩性误差。由于该方法模拟的是弱可压缩流体，当马赫数过大时，密度波动会显著影响速度场精度。其次，松弛时间必须大于二分之一，以保证正黏度。当雷诺数增大时，松弛时间趋近于二分之一，碰撞过程变慢，数值振荡容易累积并导致发散。此外，网格分辨率与时间步长之间需满足特定关系：当网格加密时，时间步长需相应减小，以保持马赫数在合理范围内。

7.3 高雷诺数失稳机理

在高雷诺数下，BGK碰撞算子的稳定性显著下降。其根本原因在于：高雷诺数对应低黏度，低黏度对应松弛时间接近二分之一，此时分布函数偏离平衡态的程度较大，线性松弛假设不再成立。此外，高雷诺数下涡结构尺度减小，需要更细的网格才能分辨，进一步加剧了计算负担。为克服这一困难，可采用多松弛时间模型，对不同物理量施加不同的松弛速率，从而提高高雷诺数下的稳定性。

8. 讨论

8.1 方法优势

格子Boltzmann方法在多孔介质流动模拟中展现出多方面优势。算法结构上，碰撞与迁移均为局部操作，天然适合并行计算，可借助图形处理器实现大规模加速。边界处理上，反弹规则与布尔掩膜使得复杂几何的描述变得异常简单，无需生成贴体网格。几何建模上，可直接导入二值图像，实现从真实岩石图像到数值模拟的无缝衔接。物理透明性上，分布函数的演化直观反映流体微观动力学，便于理解与调试。

8.2 方法局限

该方法也存在若干局限。可压缩性误差方面，格子Boltzmann方法模拟不可压缩流动时存在与马赫数平方成正比的误差，需通过减小马赫数加以抑制。高雷诺数稳定性方面，BGK算子在松弛时间接近二分之一时容易失稳，需借助多松弛时间模型或大涡模拟加以改善。单相流限制方面，本文仅考虑单相流体，未涉及多相驱替与界面动力学。二维简化方面，真实多孔介质为三维结构，二维模拟可能低估孔隙连通性与各向异性对流动的影响。

8.3 工程应用前景

尽管存在上述局限，格子Boltzmann方法在孔隙尺度流动研究中仍具有广阔的应用前景。在油气开发中，可用于预测相对渗透率曲线、模拟剩余油分布、评估提高采收率措施的微观机理。在水文地质中，可用于模拟污染物在含水层中的迁移路径与滞留时间。在化学工程中，可用于优化多孔催化剂与填充床反应器的结构设计。在生物医学中，可用于模拟血液在微血管网络中的流动与药物递送。随着计算机硬件的发展与数值方法的改进，格子Boltzmann方法有望在更大尺度、更高雷诺数与更多相态的问题中发挥更大作用。

9. 结论与展望

本文围绕二维通道内多孔介质流动问题，系统构建了基于格子Boltzmann方法的数值模拟框架。研究采用九速度离散模型与BGK碰撞算子，结合Zou–He边界条件与反弹浸入边界，依次模拟了Poiseuille流动、圆柱与方柱绕流、裂缝介质流动以及正弦孔道流动。通过与解析解和文献基准数据的对比，验证了模型在低雷诺数下的精度与可靠性。收敛性分析表明该格式具有接近二阶的空间精度；稳定性分析揭示了马赫数、松弛时间与雷诺数之间的制约关系。

未来工作可从以下方向展开：引入多松弛时间碰撞算子以提高高雷诺数稳定性；扩展至三维十九速度或二十七速度模型以捕捉真实孔隙结构；发展多相流模型以模拟油水驱替与界面动力学；结合图形处理器并行计算以加速大规模模拟；探索机器学习方法从图像中自动提取孔隙网络参数并构建升尺度模型。上述工作的推进，将进一步提升格子Boltzmann方法在多孔介质流动研究中的实用价值与科学意义。

📚第二部分——运行结果

部分代码

%%%%%%%%%%%%%%%%% INPUT 输入参数设置 %%%%%%%%%%%%%%%%%%%%%%%%%
clear % 清空工作区变量
griditer=1; % 网格迭代计数器，用于不同网格尺寸算例
% 循环遍历3套网格：ny分别取100、200、400，nx=ny，正方形网格
for ny = [100 200 400]
nx = ny; % x方向网格数等于y方向网格数
Re = 100; % 雷诺数
U = .1; % 入口最大速度
%%%%%%%%%%%%%%%% IMMERSED OBJECT 浸入式障碍物（固壁）%%%%%%%%%%%%%%
obst = zeros(nx,ny); % 障碍物标记矩阵，0=流体域，1=固体壁面
[y,x] = meshgrid(1:ny,1:nx); % 生成网格坐标矩阵
obst(:,[1,ny]) = 1; % 将左右两列标记为固体边界
immersed = find(obst);% 找出所有固体网格点的线性索引
figure;imshow(obst) % 绘图展示障碍物分布
%%%%%%%%%%%%%%%% OUTER BOUNDS 外边界与入口速度剖面 %%%%%%%%%%%%%%%%%
L = ny-2; % 流道有效长度
H = y-nx/ny;
% 抛物线型入口速度分布（泊肃叶剖面）
u = 4*U/(L^2)*(H*L-H.^2);
v = zeros(nx,ny); % 初始y方向速度置0
%%%%%%%%%%%%%%%% INSTANTIATION LBM模型初始化 %%%%%%%%%%%%%%%%%
tau = 0.5+3*U*ny/Re; % 松弛时间，由雷诺数计算
%tau = 0.5+6*U*r/Re; % 备用松弛时间公式，本算例未启用
w = [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]; % D2Q9权系数
cx=[0 1 0 -1 0 1 -1 -1 1]; % D2Q9离散速度x分量
cy=[0 0 1 0 -1 1 1 -1 -1]; % D2Q9离散速度y分量
oppdir=[1 4 5 2 3 8 9 6 7];% 反弹边界：当前方向对应的反方向编号
c=[cx;cy]'; % 离散速度向量矩阵
% u = zeros(ny,nx); % 备用：速度初始化，本代码不使用
% v = zeros(ny,nx);
%rho = ones(ny,nx);
rho = 1; % 流体密度初始值
feq=zeros(9,nx,ny); % 平衡态分布函数 9个方向 × nx × ny网格
f=zeros(9,nx,ny); % 分布函数
fnx=zeros(9,nx,ny); % 碰撞后待迁移的分布函数
% 计算D2Q9平衡态分布函数
for i=1:9
feq(i,:,:)=w(i)*rho.*(1+3*(c(i,1)*u+c(i,2)*v)+ 9/2*(c(i,1)*u+c(i,2)*v).^2 - 3/2*(u.^2+v.^2));
end;
f = feq; % 初始分布函数等于平衡态
iter=1; % 时间迭代步计数器
uerror=100; % u速度收敛误差初值
verror=100; % v速度收敛误差初值
%%%%%%%%%%%%%%%%%%%%%%%% START ITERATIONS 迭代主循环 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 迭代终止条件：迭代>1000 或 u误差<0.001 且 v误差<0.001
while(iter<1000 || uerror>.001 || verror >.001)
uprev=reshape(u,nx,ny); % 保存上一步u速度，用于收敛判断
vprev=reshape(v,nx,ny); % 保存上一步v速度
%%%%%%%%%%%%%%%%%%%%%%%% MACRO VARIABLES 宏观量（密度、速度）求解 %%%%%%%%%%%%%%%%%%%%%%%%
rho = sum(f); % 密度：所有方向分布函数求和
% 下面注释块：原始逐循环计算宏观动量，改用矩阵乘法加速
% rhou=f(1,:,:)*c(1,1);
% for i=2:9
% rhou=rhou+f(i,:,:)*c(i,1);
% end;
% rhov=f(1,:,:)*c(1,2);
% for i=2:9
% rhov=rhov+f(i,:,:)*c(i,2);
% end;
% u=squeeze(rhou./rho);
% v=squeeze(rhov./rho);
% rho=squeeze(rho);
% 矩阵运算快速计算宏观速度u,v
u=reshape((cx * reshape(f,9,nx*ny)),1,nx,ny)./rho;
v=reshape((cy * reshape(f,9,nx*ny)),1,nx,ny)./rho;
%%%%%% BOUNDARIES 边界处理
% INLET 入口边界（左边界）
H =[2:(ny-1)]-1.5;
u(:,1,2:(ny-1))=4*U/(L^2)*(H*L-H.^2); % 入口抛物线速度剖面
v(:,1,[2:(ny-1)])= 0; % 入口y向速度为0
% 入口密度重构
rho(:,1,[2:(ny-1)])=1./(1-u(:,1,[2:(ny-1)])).*(sum(f([1,3,5],1,[2:(ny-1)]))+2*sum(f([4,7,8],1,[2:(ny-1)])));
% 入口未知分布函数更新（Zou-He速度边界格式）
f(2,1,[2:(ny-1)])=f(4,1,[2:(ny-1)])+2/3*rho(:,1,[2:(ny-1)]).*u(:,1,[2:(ny-1)]);
f(6,1,[2:(ny-1)])=f(8,1,[2:(ny-1)])+1/2*(f(5,1,[2:(ny-1)])-f(3,1,[2:(ny-1)]))+ 1/2*rho(:,1,[2:(ny-1)]).*v(:,1,[2:(ny-1)]) ...
+ 1/6*rho(:,1,[2:(ny-1)]).*u(:,1,[2:(ny-1)]);
f(9,1,[2:(ny-1)])=f(7,1,[2:(ny-1)])+1/2*(f(3,1,[2:(ny-1)])-f(5,1,[2:(ny-1)]))-1/2*rho(:,1,[2:(ny-1)]).*v(:,1,[2:(ny-1)]) ...
+ 1/6*rho(:,1,[2:(ny-1)]).*u(:,1,[2:(ny-1)]);
% OUTLET 出口边界（右边界）
rho(:,nx,[2:(ny-1)])=1; % 出口密度设为1
u(:,nx,[2:(ny-1)])=-1+1./(rho(:,nx,[2:(ny-1)]))*(sum(f([1,3,5],nx,[2:(ny-1)]))+2*sum(f([2,6,9],nx,[2:(ny-1)])));
v(:,nx,[2:(ny-1)])=0; % 出口y向速度0
% 出口未知分布函数更新（Zou-He压力/出口边界）
f(4,nx,[2:(ny-1)])=f(2,nx,[2:(ny-1)])-2/3*rho(:,nx,[2:(ny-1)]).*u(:,nx,[2:(ny-1)]);
f(8,nx,[2:(ny-1)])=f(6,nx,[2:(ny-1)])+1/2*(f(3,nx,[2:(ny-1)])-f(5,nx,[2:(ny-1)]))- 1/2*rho(:,nx,[2:(ny-1)]).*v(:,nx,[2:(ny-1)]) ...
- 1/6*rho(:,nx,[2:(ny-1)]).*u(:,nx,[2:(ny-1)]);
f(7,nx,[2:(ny-1)]) = f(9,nx,[2:(ny-1)]) + 1/2*(f(5,nx,[2:(ny-1)])-f(3,nx,[2:(ny-1)]))+ 1/2*rho(:,nx,[2:(ny-1)]).*v(:,nx,[2:(ny-1)]) ...
- 1/6*rho(:,nx,[2:(ny-1)]).*u(:,nx,[2:(ny-1)]);
%%%%%%%%%%%%%%%%%%%% F EQIILIBRIUM 重新计算平衡态分布函数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
collOperator=(1-1/tau)*f+1/tau*feq;
for i=1:9
feq(i,:,:)=w(i)*rho.*(1+3*(c(i,1)*u+c(i,2)*v)+ 9/2*(c(i,1)*u+c(i,2)*v).^2 - 3/2*(u.^2+v.^2));
fnx(i,:,:)=f(i,:,:)-(1/tau).*(f(i,:,:)-feq(i,:,:)); % LBM碰撞步
end
%%%%%%%%%%%%%%%%%%%% COLLISION 碰撞+固体反弹边界 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%BOUNCE-BACK 固体壁面反弹边界条件
for i=1:9
fnx(i,immersed)=f(oppdir(i),immersed);
end
%%%%%%%%%%%%%%%%% PROPAGATION 迁移步（streaming） %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:9
f(i,:,:)=circshift(fnx(i,:,:),[0,cx(i),cy(i)]); % 分布函数沿离散速度方向迁移
end
%%%%%%%%%%%%%%%%%%% STABILITIY CRITERIA 收敛误差计算 %%%%%%%%%%%%%%%%%%%%%%%
% 计算u场Frobenius范数误差
uerror=norm((reshape(u,nx,ny)-uprev)/nx/ny,'fro');
verror=norm((reshape(u,nx,ny)-uprev)/nx/ny,'fro'); % 【注意：源码这里疑似bug，应该是vprev】
iter=iter+1
end
%%%%%%%%%%%%%%%%%% PLOTS 结果后处理绘图 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
finError(griditer)=uerror; % 保存当前网格收敛误差
griditer=griditer+1;
end
grid =[10,20,40];
% 计算空间收敛阶（网格收敛性分析）
spatial_error=-log((finError(2)-finError(3))/(finError(1)-finError(2)))/log(2)
% 生成绘图坐标
x=.5/nx:1/nx:1;
y=.5/ny:1/ny:1;
u(immersed) = nan; % 固体点置nan，绘图不显示
v(immersed) = nan;
UX=reshape(u,nx,ny);
UY=reshape(v,nx,ny);
figure;streamslice(x,y,UX',UY') % 绘制流线切片
title('Streams')
xlabel('x')
ylabel('y')
axis tight
[X,Y]=meshgrid(x,y);
[sx,sy] = meshgrid(.05,.05:.05:.95);
figure;
streamline(stream2(X,Y,UX',UY',sx,sy)) % 绘制流线
title('Streamlines')
xlabel('x')
ylabel('y')
% 绘制速度幅值云图
figure;contourf(x,y,hypot(UX',UY'))

🎉第三部分——参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。(文章内容仅供参考，具体效果以运行结果为准)

​​​​​​🌈第四部分——本文完整资源下载

资料获取，更多粉丝福利，MATLAB|Simulink|Python|数据|文档等完整资源获取

本文完整资源下载