第３６卷第３期 Ｖｏ１．３６ ＮＯ．３ 宁夏大学学报（自然科学版） Ｊｏｕｒｎａｌ ｏｆ Ｎｉｎｇｘｉａ Ｕｎｉｖｅｒｓｉｔｙ（Ｎａｔｕｒａｌ Ｓｃｉｅｎｃｅ Ｅｄｉｔｉｏｎ） ２Ｏ１５年９月 Ｓｅｐ．２０１５ 文章编号：０２５３－２３２８（２０１５）０３—０２５１—０４ 基于格子玻尔兹曼方法的自仿射粗糙裂隙渗流模拟 吴 伟， 张晓哲， 詹光亮 （宁夏回族自治区电力设计院，宁夏银川 ７５０００１） 摘 要：基于格子玻尔兹曼方法模拟研究了自仿射粗糙裂隙的渗流特征．通过连续随机叠加技术，生成了不同粗糙 指数的自仿射粗糙裂隙．为验证格子玻尔兹曼方法，将格子玻尔兹曼结果与泊肃叶流解析解对比．通过格子玻尔兹 曼方法与自仿射粗糙裂隙技术的结合，成功模拟了粗糙裂隙内部的流场分布，并分析了粗糙裂隙内部的渗流路径 与水力等效隙宽的变化．结果表明，不均匀的隙宽分布对裂隙的流场分布影响较大；尽管两组裂隙内部的Ｒｅ一致， 但其最大流速有一定差异． 关键词：格子玻尔兹曼方法；自仿射；粗糙裂隙；渗流 分类号：（中图）Ｏ３５ 文献标志码：Ａ 随着现代工业的快速发展和世界人口数量的激 常良好的并行特性及复杂边界适应性强和计算程序 编写简单等优点＿８］． 增，人们对地下水资源的需求量越来越大［１］．裂隙介 质作为地下水最主要和最重要的载体，裂隙介质中 的地下水渗流问题一直受到广泛关注 ］．通常，天然 裂隙介质结构是十分复杂的．早期，研究人员均以理 想光滑的裂隙介质为模型研究渗流问题．著名的裂 隙水流立方定律（ｃｕｂｉｃ ｌａｗ）就是假设裂隙壁为理想 光滑的，并根据大量实验而得到的． 本文通过格子玻尔兹曼方法，模拟计算了不同 粗糙程度的单裂隙渗流运移，采用连续随机叠加技 术，生成具有不同粗糙指数的自仿射粗糙裂隙．通过 理想单裂隙模型，验证了格子玻尔兹曼方法的准确 性和有效性．将自仿射粗糙裂隙与格子玻尔兹曼方 法相结合，分析了粗糙程度对渗流运移的影响． 格子玻尔兹曼方法（１ａｔｔｉｃｅ ｂｏｌｔｚｍａｎｎ ｍｅｔｈｏｄ） 是一种基于气体动理论的数值模拟方法［３ ］．格子 玻尔兹曼方法与传统计算流体力学的数值模拟方法 （ｃｏｍｐｕｔａｔｉｏｎａｌ ｆｌｕｉｄ ｄｙｎａｍｉｃｓ，ＣＦＤ）不同，格子玻 尔兹曼方法在宏观上是离散方法，而在微观上是连 １格子玻尔兹曼方法 本文利用单松弛因子（ｒ）的离散格子玻尔兹曼 方程表示裂隙介质中的流体．在微观层面中，相应计 算节点上的粒子分布函数方程为 厂 （Ｘ＋Ｐ △ ，ｔ＋△ ）一 续方法，因而通常被称为介观模拟方法．格子玻尔兹 曼方法的基本思想是认为宏观上的流体流动是大量 微观流体粒子无规则运动的结果．ＣＦＤ从微观尺度 ｆｌ（Ｘ，ｔ）一 Ｅｆｉ（Ｘ，ｔ）一ｆ （Ｘ， ）］， （１） 来研究流体运动，认为连续的流体是由大量微观粒 子构成的无序离散系统，系统中每个流体粒子都做 无规则的热运动，粒子间通过相互碰撞交换动量和 式中：ｆ （Ｘ， ）为ｔ时刻计算节点Ｘ上ｉ方向的粒 子分布函数；ｉ为格子离散速度的第ｉ个方向； 为ｉ 方向上的离散格子速度；△ｚ为计算节点的格子长 度；△ 为计算的时间步长．本文选择Ｄ２Ｑ９离散速 度模型（图１），其中，Ｐ 共有９个速度分量，即 能量，近而表现出宏观上的流体流动行为．近年来， ＣＦＤ被广泛应用到许多传统模拟方法难以胜任的 研究领域，例如：界面科学、传热机理、生物流体、材 ［ ，ｅｌ，ｅ２， ， ４， ５，ｅ６，ｅ７，ｅ８］＝＝＝ ｒ０ １ ０—１ ０ ０ １—１１ —１ １］ 料晶体生长、传质机理、计算机科学等 ］，并取得 很大成功．此外，ＣＦＤ的离散格式方程形式具有非 收稿日期：２０１４—０８—０８ Ｃｌ Ｌ０ ０ ］一】 】 １．（２） —１ —１Ｊ 作者简介：吴伟（１９８６一），男，硕士，主要从事裂隙介质水方面的研究 ２５２ 宁夏大学学报（自然科学版） 第３６卷 式中：Ｃ一△ ／Ａｔ．（１）式中，ｆ。 （Ｘ，￡）表示粒子的 直接生成，其隙宽分布符合分形理论．大量研究资 料表明，自然粗糙裂隙面是一种自仿射分形面，其 相应的高程分布符合高斯分布．二维自仿射分布 平衡分布函数，其表达式为 ． （Ｘ， ）一 ［ ＋ ＋ 一 ］，㈦ 变换分别为 满足标度不变特性，其相应的３：７和Ｙ方向的仿射 式中：ｃ 为格子固有常数，一般ｃ 一１／√３； 。为平 衡分布函数的速度； 为权重系数，对于Ｄ２Ｑ９模 ．Ｚ ．Ｚ ．Ｚ 型而言，其相应的表达式为 一 一 一 Ｏ １ ５ ｆ４／９， ２ ６ 一 １／９， ３ ７ ，●， ４ ， ８● 【１／３６。 ， 、，４ 图１ Ｄ２Ｑ９离散速度模型 流体的宏观物理参数可通过分析计算节点上微 观粒子分布函数而得到，相应的表达式如下： Ｐ一∑厂 ， （５） ０ Ｐ—ｃ ２ ， （６） ｕ—Ｃ：（ｒ一０．５）Ａｔ． （７） 对于格子玻尔兹曼方法，其边界条件种类十分 多样．常见的有：标准反弹边界、对称边界、充分发展 边界和周期性边界等．在本文所有模拟中，计算区域 中的节点只有两种，即固体节点和流体节点．对于固 体节点，其边界条件均为标准反弹边界．计算中，通 过在每一个流体节点施加体积力Ｇ，使流体发生流 动．在计算区域的进口和出口处，采用周期性边界， 这样可使整个区域的质量保持守恒． ２ 白仿射粗糙裂隙的构建 据现有的大量资料证明，天然开裂的粗糙裂 隙表面具有分形特征．这样，便可以在计算机上实 现对自然开裂的裂隙面粗糙度及隙宽分布的构 建．Ａｕｒａｄｏｕ等 通过有机玻璃制作了自仿射粗糙 单裂隙，并在其裂隙内研究不可混溶两相流体的 运移问题，该自仿射粗糙单裂隙通过统计学特征 ｚ 一 ， （８） 一Ｙ ， （９） 式中： 被称为粗糙指数或Ｈ指数（Ｈｕｓｒｔ Ｅｘｐｏｎｅｎｔ）．本文应用的二维自仿射线生成技术是 由Ｖｏｓｓ＿ｌ叩于１９８５年提出的，其生成技术的基本思 路是依据（８）和（９）式，通过而且符合正态高斯 分布的变量 ，对选定点逐个随机生成高程，并通过 不断插值，获得两点之间的新高程．二维自仿射线的 生成过程如下： 首先，在平均数（Ｍｅａｎ）为０、方差（Ｖａｒｉａｎｃｅ）为 １符合高斯分布的集合 中，通过抽取一点作为 一 １／２点上的高程值．高程值即为ｙ方向上的值． 然后，通过ｚ一０和３２—１／２的点上的高程插值 获得ｚ一１／４点上的高程，同理，ｚ一３／４的点上的高 程通过ｚ＝１／２和ｚ一１的点上的高程插值获得．这 样，对于新生成的高程随机序列，其方差为 斋‘ Ｖｏｓｓ 通过功率谱分析证明，该算法生成的曲 线符合自仿射特性，即在３２和Ｙ方向上，符合（８）和 （９）式．图２显示了在Ｈ指数为０．８１时，二维白仿 射曲线的示意图． 图２粗糙指数为０．８１时的二维自仿射曲线 ３分析与讨论 ３．１模拟验证 为验证格子玻尔兹曼方法的准确性，本文采用 理想光滑的单裂隙作为验证算例．对于理想单裂隙 模型，假设水流流动状态为层流并考虑实际流体的 黏滞力，这样的流体流动通常被称为泊肃叶流 （Ｐｏｉｓｅｕｉｌｌｅ’Ｓ ｆｌｏｗ）．在水平方向上的速度分布，即 流体速度在隙度方向上的分布，可由（２）式得 第３期 吴 伟等：基于格子玻尔兹曼方法的自仿射粗糙裂隙渗流模拟 ２５３ （ ）一一 ｂ２ ｆ１一 １， （１１） 部流场达到稳定．本文所有单位均为格子玻尔兹曼单 ６“ｄＺ＼０， 式中：“为水流速度；ｂ为理想光滑裂隙的隙宽；５Ｃ为 水平方向坐标； 为垂直方向（隙宽方向）坐标； 为 流体动力黏滞系数．本文取理想光滑单裂隙的大小 为２５６×６０ ｌｕ，裂隙壁采用标准反弹边界条件，裂隙 的出口和进口均采用周期边界，通过施加体积力 Ｇ一０．０００ ０２，使裂隙内流体流动．图３为理想光滑 单裂隙内的流场分布，其中：红色代表裂隙壁；蓝色 代表裂隙内被水充填的区域；灰白色的箭头为流场． 图３理想光滑裂隙内部流场 从图３中可以看出，理想光滑单裂隙内的流场 分布并不一致，隙宽中间区域的流速较大，靠近裂隙 壁处的流速较小．这一结果与（１１）式一致．截取光滑 裂隙隙宽方向上的速度分布，与（１１）式进行比较，结 果见图４． 图４隙宽方向上流速解析解与ＬＢＭ结果 ３．２粗糙裂隙渗流特征 本文通过４条不同粗糙指数的自仿射曲线组成２ 组粗糙裂隙．图５和６为利用格子玻尔兹曼方法计算 的粗糙裂隙稳定流场．对于图５和６，其计算区域的大 小为６０Ｘ２６０ ｌｕ ，裂隙内的水流方向为从左至右．图５ 和６中的边界条件设置相同，对于粗糙裂隙壁，采用标 准反弹边界，粗糙裂隙的进口和出口，均采用周期边界 条件．通过施加体积力Ｇ使粗糙裂隙内的流体流动．分 别调整体积力的大小，使图５和６两组粗糙裂隙内的 Ｒｅ一１５．因此，图５和６中的流体流动均属层流．一般 情况，利用格子玻尔兹曼方法时，粗糙裂隙内的流场需 要迭代多次才能到达稳定．本文假设粗糙裂隙内部的 平均流速在５ ０００ ｔｓ内的误差小于０．１ ９／６时，粗糙裂隙内 位，例如：长度为１ｕ，时间为ｔｓ等． 图５中粗糙裂隙上方自仿射曲线的粗糙指数为 ０．８２，下方的粗糙指数为０．７３；图６中粗糙裂隙上 方自仿射曲线的粗糙指数为０．７８，下方的粗糙指数 为０．８５． 对比图５，６与图３的流场差异，可以发现，对于 粗糙的裂隙壁，其内部流场分布并不均匀．在靠近粗 糙裂隙壁处，流体流速较小，而在粗糙裂隙中部，流 体流速较大，这符合实际情况．尽管２组粗糙裂隙的 Ｒｅ－－１５，但是由于粗糙裂隙壁造成的隙宽分布不均 匀，图５和６的内部流速不一致．图５和６中裂隙内 部最大流速分别为０．００１ ３９８和０．０１８ ５ ｌｕ／ｔｓ．从 图５和６可以看出，在隙宽较小的位置，图６中的流 体流速明显大于图５中的流体流速． 图５ 平均隙宽为３５ ｌｕ的粗糙裂隙流场分布 图６ 平均隙宽为２１ ｌｕ的粗糙裂隙流场分布 ２组粗糙裂隙的平均隙宽分别为３５ ｌｕ（图５）和 ２１ ｌｕ（图６）．平均隙宽是粗糙裂隙在水平方向上各 点隙宽值的几何平均数．然而，在实际问题中，平均 隙宽较难获得，需要通过立方定律（ｃｕｂｉｃ ｌａｗ）反算 求得水力等效隙宽．对于理想光滑裂隙而言，求得的 水力等效隙宽等于平均隙宽，但对于粗糙裂隙而言， 其水力等效隙宽与平均隙宽通常存在一定的差别． 立方定律的表达式如下： ３ ｑ一曲一 ， （１２） 上 ｗ 式中：ｑ为裂隙内的流量； 为水流平均速度；ｇ为当 地重力加速度； 为水的运动粘滞系数； 为与裂隙 平面平得的水力坡降．需要指出的是，（１２）式中的水 力坡降 ，可以用格子玻尔兹曼方法中的体积力Ｇ 代替．由于本文所有单位均为格子玻尔兹曼单位，因 而需要对当地的重力加速度ｇ进行量纲转换． 通过多次改变体积力，并将获得的多组隙宽进 行几何平均，最后获得水力等效隙宽．利用（１２）式对 图５和６求得的水力等效隙宽分别为３０．２ ｌｕ和 ２５４ 宁夏大学学报（自然科学版） 第３６卷 １４ ｌｕ．可以发现，两者的水力等效隙宽均小于其平 均隙宽．并且，图５的平均隙宽与水力等效隙宽的差 ｗｅｔｔａｂｉｌｉｔｙ ｏｎ ｉｎｔｅｒｒａｃｉａｌ ａｒｅａ ｄｕｒｉｎｇ ｉｍｍｉｓｃｉｂｌｅ ｌｉｑｕｉｄ ｉｎｖａｓｉｏｎ ｉｎｔｏ ａ ３Ｄ ｓｅｌｆ—ａｆｆｉｎｅ ｒｏｕｇｈ ｆｒａｃｔｕｒｅ：ｌａｔｔｉｃｅ 值小于图６的差值．这反映了图６的裂隙粗糙程度 要大于图５的裂隙粗糙程度． ｂｏｌｔｚｍａｎｎ ｓｉｍｕｌａｔｉｏｎｓ［Ｊ］． Ａｄｖａｎｃｅｓ ｉｎ Ｗａｔｅｒ Ｒｅｓｏｕｒｃｅｓ，２０１３（６１）：１—１１ ［３］郭照立，郑楚光．格子Ｂｏｌｔｚｍａｎｎ方法的原理及应用 ４ 结论 １）格子玻尔兹曼方法是研究粗糙裂隙内部的流 场的有效方法，不仅可以模拟真实的粗糙裂隙内部 流场分布，而且具有处理复杂边界的能力． ［Ｍ］．北京：科学出版社，２００９ ［４］何雅玲，王勇，李庆．格子Ｂｏｌｔｚｍａｎｎ方法的理论及应 用［Ｍ］．北京：科学出版社，２００９ Ｅ５２ 吴伟，孙东科，戴挺，等．枝晶生长和气泡形成的数值模 拟［Ｊ］．物理学报，２０１２，６１（１５）：１－９ ［６］ 孙东科，项楠，陈科，等．格子玻尔兹曼方法模拟弯流道中 粒子的惯性迁移行为［Ｊ］．物理学报，２０１３，６２（２）：１－９ ２）不均匀的隙宽分布对裂隙的流场分布影响较 大．虽然两组裂隙内部的Ｒｅ一致，但其最大流速有 一定差异．这种差异受到了裂隙的粗糙程度影响，粗 ３）通过立方定律反算的水力等效隙宽小于裂隙 ［７３孙东科．格子Ｂｏｌｔｚｍａｎｎ方法模拟合金对流枝晶的生 长［ｈｉ．南京：东南大学，２０１０ ［８］ＳＵＣＣＩ Ｓ．Ｔｈｅ ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ ｅｑｕａｔｉｏｎ：ｆｏｒ ｆｌｕｉｄ ｄｙｎａｍｉｃｓ ａｎｄ ｂｅｙｏｎｄ［Ｍ］．Ｏｘｆｏｒｄ：Ｏｘｆｏｒｄ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ，２００１． 糙程度越大，速流差异越大． 的平均隙宽，两者之间的差值可以反映裂隙的粗糙 程度差异． ［９］ＡＵＲＡＤＯＵ Ｈ，ＨＵＬＩＮ Ｊ Ｐ，ＲＯＵＸ Ｓ．Ｅｘｐｅｒｉｍｅｎｔａｌ 参考文献 ［１］ 窦智，周志芳，李兆峰．多孔介质油水两相ｋ～ｓ～Ｐ关 系数学模型的实验研究ＥＪ］．水科学进展，２０１２，２３（２）： ２０６—２１３． ｓｔｕｄｙ ｏｆ ｍｉｓｃｉｂｌｅ ｄｉｓｐｌａｃｅｍｅｎｔ ｆｒｏｎｔｓ ｉｎ ｒｏｕｇｈ ｓｅｌｆ ａｆｆｉｎｅ ｆｒａｃｔｕｒｅｓ［Ｊ］．Ｐｈｙｓｉｃａｌ Ｒｅｖｉｅｗ Ｅ，２００１，６３ （６）：０６６３０６． ［１ｏ］ ＶＯＳＳ Ｒ Ｆ．Ｒａｎｄｏｍ ｆｒａｃｔａｌ ｆｏｒｇｅｒｉｅｓ［Ｊ］．Ｆｕｎｄａｍｅｎｔａｌ Ａｌｇｏｒｉｔｈｍｓ ｆｏｒ Ｃｏｍｐｕｔｅｒ Ｇｒａｐｈｉｃｓ，１９９１（１７）：８０５—８３５ ［２］ＤＯＵ Ｚｈｉ，ＺＨＯＵ Ｚｈｉｆａｎｇ，ＳＬＥＥＰ Ｂ Ｅ．Ｉｎｆｌｕｅｎｃｅ ｏｆ Ｓｅｌｆ－ａｆｆｉｎｅ Ｒｏｕｇｈ Ｆｒａｃｔｕｒｅ Ｆｌｏｗ Ｓｉｍｕｌａｔｉｏｎ Ｂａｓｅｄ ｏｎ ｔｈｅ Ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ Ｍｅｔｈｏｄ Ｗｕ Ｗｅｉ，Ｚｈａｎｇ Ｘｉａｏｚｈｅ，Ｚｈａｎ Ｇｕａｎｇｌｉａｎｇ （Ｎｉｎｇｘｉａ Ｅｌｅｃｔｒｉｃ Ｐｏｗｅｒ Ｄｅｓｉｇｎ Ｉｎｓｔｉｔｕｔｅ，Ｙｉｎｃｈｕａｎ ７５０００１，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ｔｈｅ ｓｅｌｆ－ａｆｆｉｎｅ ｒｏｕｇｈ ｆｒａｃｔｕｒｅ ｓｅｅｐａｇｅ ｆｅａｔｕｒｅ ｉｓ ｓｔｕｄｉｅｄ ｂｙ ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ ｍｅｔｈｏｄ．Ａ ｓｅｌｆ－ ａｆｆｉｎｅ ｒｏｕｇｈ ｆｒａｃｔｕｒｅ ｗｉｔｈ ｄｉｆｆｅｒｅｎｔ ｉｎｄｉｃｅｓ ｉｓ ｇｅｎｅｒａｔｅｄ ｔｈｒｏｕｇｈ ｃｏｎｔｉｎｕｏｕｓ ｒａｎｄｏｍ ｓｔａｃｋｉｎｇ ｔｅｃｈｎｉｑｕｅ．Ｔｏ ｖｅｒｉｆｙ ｔｈｅ ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ ｍｅｔｈｏｄ。ｔｈｅ ｒｅｓｕｌｔｓ ｏｆ ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ ａｎｄ ａｎａｌｙｔｉｃａｌ ｓｏｌｕｔｉｏｎ ｏｆ Ｐｏｉｓｅｕｉｌｌｅ ｆｌｏｗ ａｒｅ ｃｏｎｔｒａｓｔｅｄ．Ｂｙ ｃｏｍｂｉｎｉｎｇ ｓｅｌｆ－ａｆｆｉｎｅ ｒｏｕｇｈ ｆｒａｃｔｕｒｅ ｔｅｃｈｎｉｑｕｅ ａｎｄ ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ ｍｅｔｈｏｄ，ｔｈｅ ｆｌｏｗ ｆｉｅｌｄ ｄｉｓｔｒｉｂｕｔｉｏｎ ｉｎｓｉｄｅ ｔｈｅ ｒｏｕｇｈ ｆｒａｃｔｕｒｅ ｉｓ ｓｉｍｕｌａｔｅｄ ｓｕｃｃｅｓｓｆｕｌｌｙ，ａｎｄ ｔｈｅ ｒｏｕｇｈ ｆｒａｃｔｕｒｅ ｓｅｅｐａｇｅ ｐａｔｈ ｉｎｓｉｄｅ ｔｈｅ ｈｙｄｒａｕｌｉｃ ｅｑｕｉｖａｌｅｎｔ ｇａｐ ｗｉｄｔｈ ｃｈａｎｇｅｓ ｉｓ ａｎａｌｙｚｅｄ．Ｔｈｅ ｒｅｓｕｌｔｓ ｓｈｏｗ ｔｈａｔ ｔｈｅ ｕｎｅｖｅｎ ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ ｇａｐ ｗｉｄｔｈ ｈａｓ ｇｒｅａｔｅｒ ｉｍｐａｃｔ ｏｎ ｔｈｅ ｆｒａｃｔｕｒｅ ｆｌｏｗ ｆｉｅｌｄ ｄｉｓｔｒｉｂｕｔｉｏｎ．Ｄｅｓｐｉｔｅ Ｒｅ ｎｕｍｂｅｒ ｉｎｓｉｄｅ ｔｈｅ ｒｏｕｇｈ ｆｒａｃｔｕｒｅ ｉＳ ｃｏｎｓｉｓｔｅｎｔ。ｔｈｅｒｅ ａｒｅ ｓｏｍｅ ｄｉｆｆｅｒｅｎｃｅｓ ｉｎ ｔｈｅ ｍａｘｉｍｕｍ ｆｌｏｗ ｒａｔｅ ｂｅｔｗｅｅｎ ｔＷＯ ｆｒａｃｔｕｒｅｓ． Ｋｅｙ ｗｏｒｄｓ：ｌａｔｔｉｃｅ Ｂｏｌｔｚｍａｎｎ ｍｅｔｈｏｄ；ｓｅｌｆ—ａｆｆｉｎｅ；ｒｏｕｇｈ ｆｒａｃｔｕｒｅ；ｓｅｅｐａｇｅ （责任编辑、校对张刚）
