欢迎来到本博客❤️❤️博主优势博客内容尽量做到思维缜密逻辑清晰为了方便读者。完整资源、论文复现、期刊合作、论文辅导及科研仿真定制事宜点击本文完整资源下载⛳️座右铭行百里者半于九十。⛳️赠与读者‍做科研涉及到一个深在的思想系统需要科研者逻辑缜密踏实认真但是不能只是努力很多时候借力比努力更重要然后还要有仰望星空的创新点和启发点。建议读者按目录次序逐一浏览免得骤然跌入幽暗的迷宫找不到来时的路它不足为你揭示全部问题的答案但若能解答你胸中升起的一朵朵疑云也未尝不会酿成晚霞斑斓的别一番景致万一它给你带来了一场精神世界的苦雨那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。或许雨过云收神驰的天地更清朗.......第一部分——内容介绍基于格子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 % 清空工作区变量griditer1; % 网格迭代计数器用于不同网格尺寸算例% 循环遍历3套网格ny分别取100、200、400nxny正方形网格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.53*U*ny/Re; % 松弛时间由雷诺数计算%tau 0.56*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; % 流体密度初始值feqzeros(9,nx,ny); % 平衡态分布函数 9个方向 × nx × ny网格fzeros(9,nx,ny); % 分布函数fnxzeros(9,nx,ny); % 碰撞后待迁移的分布函数% 计算D2Q9平衡态分布函数for i1:9feq(i,:,:)w(i)*rho.*(13*(c(i,1)*uc(i,2)*v) 9/2*(c(i,1)*uc(i,2)*v).^2 - 3/2*(u.^2v.^2));end;f feq; % 初始分布函数等于平衡态iter1; % 时间迭代步计数器uerror100; % u速度收敛误差初值verror100; % v速度收敛误差初值%%%%%%%%%%%%%%%%%%%%%%%% START ITERATIONS 迭代主循环 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 迭代终止条件迭代1000 或 u误差0.001 且 v误差0.001while(iter1000 || uerror.001 || verror .001)uprevreshape(u,nx,ny); % 保存上一步u速度用于收敛判断vprevreshape(v,nx,ny); % 保存上一步v速度%%%%%%%%%%%%%%%%%%%%%%%% MACRO VARIABLES 宏观量密度、速度求解 %%%%%%%%%%%%%%%%%%%%%%%%rho sum(f); % 密度所有方向分布函数求和% 下面注释块原始逐循环计算宏观动量改用矩阵乘法加速% rhouf(1,:,:)*c(1,1);% for i2:9% rhourhouf(i,:,:)*c(i,1);% end;% rhovf(1,:,:)*c(1,2);% for i2:9% rhovrhovf(i,:,:)*c(i,2);% end;% usqueeze(rhou./rho);% vsqueeze(rhov./rho);% rhosqueeze(rho);% 矩阵运算快速计算宏观速度u,vureshape((cx * reshape(f,9,nx*ny)),1,nx,ny)./rho;vreshape((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; % 出口密度设为1u(:,nx,[2:(ny-1)])-11./(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)*f1/tau*feq;for i1:9feq(i,:,:)w(i)*rho.*(13*(c(i,1)*uc(i,2)*v) 9/2*(c(i,1)*uc(i,2)*v).^2 - 3/2*(u.^2v.^2));fnx(i,:,:)f(i,:,:)-(1/tau).*(f(i,:,:)-feq(i,:,:)); % LBM碰撞步end%%%%%%%%%%%%%%%%%%%% COLLISION 碰撞固体反弹边界 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%BOUNCE-BACK 固体壁面反弹边界条件for i1:9fnx(i,immersed)f(oppdir(i),immersed);end%%%%%%%%%%%%%%%%% PROPAGATION 迁移步streaming %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for i1:9f(i,:,:)circshift(fnx(i,:,:),[0,cx(i),cy(i)]); % 分布函数沿离散速度方向迁移end%%%%%%%%%%%%%%%%%%% STABILITIY CRITERIA 收敛误差计算 %%%%%%%%%%%%%%%%%%%%%%%% 计算u场Frobenius范数误差uerrornorm((reshape(u,nx,ny)-uprev)/nx/ny,fro);verrornorm((reshape(u,nx,ny)-uprev)/nx/ny,fro); % 【注意源码这里疑似bug应该是vprev】iteriter1end%%%%%%%%%%%%%%%%%% PLOTS 结果后处理绘图 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%finError(griditer)uerror; % 保存当前网格收敛误差griditergriditer1;endgrid [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;UXreshape(u,nx,ny);UYreshape(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|数据|文档等完整资源获取本文完整资源下载