田宇左| Jun 姚| 雷 张| 海 孙| 文欣 杨
中国石油大学（华东）深油深气国家重点实验室，青岛 266580，中华人民共和国
**摘要**
地下氢储存（UHS）的成功在很大程度上取决于氢的注入能力，而这与多种孔隙尺度因素密切相关。由于对影响因素的研究不足，这些因素的影响及其相互作用仍不确定。本研究提出了一种伪势多组分多相格子玻尔兹曼模型来模拟含水层中的UHS，探讨了毛细数（Ca）和润湿性（θ）、与盐度相关的界面张力（IFT）变化以及由微生物活动引起的异质润湿性的协同效应。结果揭示了在不同Ca和θ组合下的流动模式转变，这些转变显著影响了储存能力。位移效率仅对IFT有轻微依赖性。微生物引起的异质润湿性具有双重效应：它可以在小孔隙中改善注入效果，但可能导致大孔隙中的优先流动，从而对最大氢饱和度产生负面影响。本研究表明，应优先考虑修改储层润湿性和优化注入速率，而形成水盐度的影响较小，不应忽视微生物引起的润湿性变化。
**引言**
由于化石燃料消耗的持续增加和环境问题的恶化，可再生二次能源，特别是氢能，受到了广泛关注[1,2]。目前，氢能的利用涉及如何大规模储存以满足长期需求。地下氢储存（UHS）提供了一种大规模、高效且安全的氢储存方式。通过利用地质构造（如盐穴、枯竭的油藏和含水层）[3],[4],[5]，UHS项目可以将氢以高压气体形式储存[6],[7],[8]。在这些地质构造中，含水层具有显著的气体储存潜力，并且在全球范围内广泛分布，被认为是UHS的最佳选择之一[3,9,10]。
含水层由含水砂岩层组成，这是一种多孔且渗透性适中的岩石构造，孔隙空间被淡水或咸水占据。当氢被注入含水层时，气体会取代水并占据部分孔隙空间，迫使水迁移。在提取氢的过程中，被排出的水会反向流动[11]，同时产生水产量。部分氢会被含水层中的水相捕获而无法回收[12,13]。此外，在储存过程中，部分氢会通过与细菌和矿物的化学反应而永久损失[14],[15],[16]。上述过程中，UHS的效果在很大程度上取决于地质构造的注入能力和储存能力，这些能力由氢-水系统的孔隙尺度流动机制决定，而这种机制非常复杂，因为氢-水的不相溶流动受到流体饱和度、粘度、界面张力、润湿性和孔隙结构等因素的影响[17],[18],[19]。因此，为了分析含水层中UHS的性能并优化氢储存效率同时最小化氢损失，研究氢-水系统的孔隙尺度多相流动至关重要。
目前，一些研究已经探讨了氢-水在孔隙尺度上的流动机制。这些研究通常关注氢-水置换过程，旨在模拟含水层中的氢注入情景[18],[20],[21],[22],[23],[24],[25],[26]。结果表明，氢通常表现为非润湿相[27],[28],[29],[30]，需要在注入过程中克服毛细压力。这种毛细效应会导致指状侵入模式，导致孔隙填充不完全和氢饱和度低[18,20,23]。当氢进入孔隙时，润湿性效应会导致氢流中断并被困在孔隙中，这一过程称为“断裂”[31,32]。这种机制会不利地影响置换过程中的氢簇传输[22,25,31]。为了提高氢注入效率，先前的研究探讨了一系列控制因素，包括位移速率、注入压力、润湿性变化、温度和孔隙结构[20,21,25,33]。主要结论表明，提高氢注入速率（与毛细数相关）和压力可以适度提高氢的注入能力[20,21]。此外，降低润湿性可以减少氢断裂事件的频率，促进连续气相的形成并提高储存能力[25]。尽管一些研究已经阐明了单个因素对多相流动的影响，但由于忽略了多个因素之间的相互作用，结果仍不确定。因此，观察到的氢饱和度随毛细数或其他因素的变化可能与传统发现不符，从而影响氢储存效率优化策略的可靠性。例如，Bagheri等人[20]观察到在高注入速率下氢饱和度意外降低，Dehury等人[21]也发现，在低毛细数条件下排水过程中的氢饱和度更高。另一方面，其他多相流领域的研究人员发现了毛细数（Ca）、粘度比（M）和润湿性（θ）在控制孔隙尺度多相流动动态方面的协同机制。研究表明，Ca和M主要决定了是发生指状流动还是均匀置换，它们的相对大小是区分流动模式的主要标准[34,35]。虽然Ca-M控制着指状流动，但Ca-θ的相互作用主要决定了孔隙空间内的流体侵入形态。研究发现，提高润湿性可以促进流体在孔隙中的协同填充，但高Ca会削弱这种效果[36,37]。通过对二氧化碳[38,39]和其他流体（如油[40]）的置换模拟，验证了Ca-θ耦合对孔隙尺度侵入模式的协同效应。因此，鉴于氢与二氧化碳或油之间的物理性质差异较大，并考虑到氢-水系统固有的高粘度比，研究Ca和θ之间的协同机制作为控制孔隙尺度多相流动的主要因素至关重要。此外，一些与UHS相关的因素尚未得到充分探索。例如，尽管现有研究已经认识到微生物引起的异质润湿性的存在，但这一关键因素尚未被纳入孔隙尺度流动研究中以量化其对异质储层中氢注入效率的影响。其他因素，如氢-水界面张力（IFT），仅在吸水过程中被研究以评估其对氢回收效率的影响[41]，而IFT对氢注入能力的影响尚未完全探索。
在上述氢-水流动研究中，研究方法可以分为实验方法和数值模拟。实验方法如岩心灌水和微流控技术在物理真实性和直接观察多相流动现象方面具有显著优势[18,22,23]。特别是与先进的可视化技术结合使用时，这些方法可以为复杂的多相流动过程提供有价值的见解。数值模拟作为一种补充方法，通过实现流动过程的完全参数化建模，不仅减少了对高分辨率测量技术的依赖，还可以系统地分离控制变量，揭示单个参数对多相流动的影响机制。更重要的是，在极端储层条件下，模拟的优势更加明显。尽管先进设备（如高压微流控技术）使得在实验中再现储层条件成为可能[42]，但高压、高温和高盐度环境（如UHS含水层）对实验硬件和成像技术提出了更高的要求。相比之下，数值模拟在这些条件下具有更大的操作灵活性和成本效益，使得在广泛理论可行条件下进行流动研究成为可能。
鉴于这些优势，本研究采用数值模拟作为主要研究工具，主要包括直接数值模拟（DNS）[43],[44],[45]、格子玻尔兹曼方法（LBM）[46,47]和孔隙网络模型（PNM）[48]。与其他方法相比，LBM在界面动力学精度、复杂孔隙结构的处理以及本质上的并行化算法框架方面具有优势[49,50]。更重要的是，对于解决氢-水系统等多相流动系统，LBM采用了一种独特的流体粒子视角，大大简化了实现复杂流体动力学的算法，特别是那些描述多相流动中分子间力的算法，并捕捉润湿性效应。这种显式解决方案避免了传统数值方法中求解多个偏微分方程时出现的离散化和模型耦合问题，从而促进了模型求解。由于这些优势，LBM已广泛应用于提高石油采收率、地下水污染、CO2封存等研究多相流动在多孔介质中的应用[51],[52],[53],[54]，同样适用于描述多孔介质中的氢-水系统多相流动[55,56]。然而，LBM在UHS研究中的应用尚未广泛[55],[56],[57]，并且关于此类应用的详细参数设置的多相流动模型报告仍然有限。根本挑战在于LBM在模拟具有大密度比和粘度比的氢-水系统流动时表现出数值不稳定性[55]。当模拟具有瞬态饱和度的UHS的置换或吸水过程时，由于边界处不同流体相的流入和流出，数值不稳定性问题更加明显[58,59]。根据已发布的单组分伪势模型，氢和水被视为使用统一状态方程（EOS）的两种相。虽然这种方法便于实现并且可以满足大密度比的要求，但也可能导致氢和水之间的非物理相变，从而影响准确性。此外，两种组分的物理性质调整依赖于相同的EOS，使得难以独立调整流体性质并难以保证热力学一致性。另一种多相流动模型——颜色梯度模型，可以有效抑制这种非物理相变，但要求氢和水的密度相似，这与流体的实际物理性质不符。因此，为了实现氢和水性质的独立调整并获得更真实的模拟结果，应开发多组分模型来描述它们的各自流动。然而，这种多组分方法引入了新的挑战，包括放大的虚假速度、热力学不一致性和计算复杂性。这些因素构成了模拟不相溶氢-水流动时必须克服的几个关键计算障碍。为了克服上述限制，本研究旨在开发一个多组分多相（MCMP）LBM框架，以便在UHS研究中更可靠地使用LBM。在这个提出的模型中，氢和水分别受各自的EOS控制，从而能够模拟复杂多孔介质结构中具有大密度比和粘度比的不相溶两相流动。
本研究首先构建并验证了一个MCMP-LBM模型，然后模拟了以下UHS情景：在第一种情景中，研究了不同润湿性和毛细数协同效应下的位移。另外两种情景分别研究了由于盐度引起的界面张力变化以及由微生物生长引起的异质润湿性。
**代码片段**
**MRT-LBM方案的伪势模型**
为了描述不相溶的MCMP流动系统，我们在D2Q9格子上应用了带有MRT方案的伪势模型，如下所示[60],[61],[62]：
$$f_{\alpha,k}(x+e_{\alpha}\delta t,t+\delta t) = f_{\alpha,k}(x,t) - (M-1)\Lambda M)\alpha\beta(f_{\beta,k}-f_{\beta,k}e_{eq}+\delta tF_{\alpha,k}'$$
其中$f_{\alpha}(x, t)$表示位置$x$和时间$t$沿离散方向$\alpha$的密度分布函数。D2Q9中的离散速度$e_{\alpha}$表示为：
$$e_{\alpha} = \begin{cases}
(0,0) & \alpha=0,
c(\cos[(\alpha-1)\pi/2],\sin[(\alpha-1)\pi/2]) & \alpha=1,2,3,4,2,
c(\cos[(2\alpha-1)\pi/4],\sin[(2\alpha-1)\pi/4]) & \alpha=5,6,7,8,
\end{cases}$$
其中$c_s^2 = c^2/3$，$c = \Delta x/\delta t$。
**多孔介质模型**
用于模拟的异质多孔介质如图2(a)所示，白色部分代表孔隙空间，黑色部分表示固体骨架。我们主要基于从Bentheimer砂岩[82],[83],[84]中提取的孔隙几何特性生成了2D连通多孔域，这种岩石类型广泛用于UHS置换实验[85],[86],[87]。模拟域的尺寸为802 lu × 799 lu。
**毛细数和润湿性的协同效应**
基于第3.2节中展示的模拟方案，分析了不同的Ca和$\theta$组合。图5展示了稳态下的流体分布和最大气体饱和度$S_g^{max}$，揭示了不同的氢指状模式。在缓慢位移条件下，孔隙内的流体填充受局部毛细力的控制。在任何给定的孔隙空间区域内，都存在一个临界毛细压力Pc，它决定了一个流体相是否会取代另一个流体相。在这项研究中，应用了MCMP LBM方法来模拟氢气注入含水层中的非均质多孔介质的过程，同时研究了毛细数（Ca）和润湿性（θ）的协同效应，以及由盐度和微生物活动引起的非均质润湿性所导致的界面张力（IFT）。主要结论总结如下：(1) 在不同的Ca和θ条件下，氢气在置换过程中表现出多种流动模式，这些模式会影响氢气的扩散过程。
作者贡献声明：
左天宇：撰写原始稿件、可视化处理、软件开发、方法论设计。
姚军：撰写与编辑、监督工作、资金筹集、概念构思。
张雷：撰写与编辑、监督工作、软件开发、方法论设计。
孙海：撰写与编辑、监督工作、软件开发。
杨文欣：撰写与编辑、验证结果、方法论设计。
利益冲突声明：
作者声明他们没有已知的财务利益冲突或个人关系，这些因素可能会影响本文所述的研究结果。
致谢：
本研究得到了国家自然科学基金（项目编号：U24B6001、52034010、42090024）和山东省自然科学基金（项目编号：ZR2022JQ23）的支持。作者感谢数字多孔介质门户网站（DPMP）提供的开源核心数据[100,101]，包括参考的一系列数据库[[82]、[83]、[84]]。
