今天,特邀博主雷庆华博士和我们一起讨论了一种新的模拟裂隙介质中全耦合流体力学过程的方法。
了解裂隙地质介质中固体变形与流体流动之间的耦合关系,对于解决地球科学和岩土工程中的许多核心问题,例如地下挖掘、油气开采、碳封存、地热生产和废物处理,具有重要的意义。本文描述了一种基于 COMSOL Multiphysics® 软件的裂隙介质中全耦合流体力学过程建模的新方法。
一般来说,模拟裂隙介质中的耦合流体力学过程存在两个主要挑战。一个是内含大量的天然裂隙的地质介质的不连续性表征,这些裂隙普遍存在许多不同的长度尺度,并经常主导系统的整体行为(参考文献 2)。另一个是流体力学耦合机制的计算,包括直接耦合(即固体和流体场之间的相互作用)和间接耦合(即岩石/裂隙性质的改变)。
在过去的几年中,已经开发了大量旨在应对这些挑战的商业软件包和开源研究代码。然而,其中大多数必须使用不同的求解器来计算流体和固体方程,因此必须通过额外的处理步骤来实现耦合,这既不方便也不高效。此外,大多数现有代码无法真正同时捕获直接和间接耦合,因此通常必须进行假设或简化。
使用 COMSOL Multiphysics 是因为它具有卓越的功能:
下面,我们将阐述在 COMSOL Multiphysics 中建立裂隙介质全耦合流体力学数值模型的步骤,并给出一些仿真实例。
在 COMSOL Multiphysics 中进行数值仿真涉及三个主要步骤。
首先,可以使用 AutoCAD® 或 Rhinoceros® 等 CAD 软件构建几何上表示为线/折线的离散裂隙网络。然后将几何数据导出为 DXFTM 文件,这些文件可以直接导入 COMSOL Multiphysics。这一步骤也可以在 MATLAB® 中完成,以按照规定的概率分布生成合成裂隙网络,并将其导出到 DXFTM。
提示:您还可以使用离散裂隙网络插件直接在 COMSOL Multiphysics 内部的现有几何结构中创建随机裂隙,如裂缝性储层的 3D 示例模型所述。
导入几何图形后,我们使用三角形有限元的非结构化网格(通过 Delaunay 细分)对域进行离散化,其中天然裂隙由嵌入相邻有限元之间的联合单元表示(图1)。
[画像:流体力学模型的离散化网格,其中标记了天然裂缝、岩石基质、节理元素、元素节点和三角形有限元]
图1。该模型采用三角形有限元的非结构化网格离散化,其中天然裂隙由嵌入相邻有限元之间的联合单元表示。
我们使用 COMSOL Multiphysics 中的固体力学 和达西定律 接口对裂隙介质中的流体力学过程进行了模拟。我们激活 多孔弹性 接口以实现固体和流体方程之间的直接耦合。我们定义了岩石基质和裂隙的材料特性和本构方程。一些岩石/裂隙特性,例如孔隙率、储水和渗透率,被定义为局部应力/压力状态的函数,以实现间接耦合。我们还定义了力学和水力边界条件。
我们在两个连续的阶段运行模型。在第一阶段,系统在给定的原位 应力和压力条件下达到初始平衡(通过斜坡加载)。然后,在第二阶段,我们模拟系统对流体注入或地下开挖等工程活动的响应。
我们应用该模型来模拟受流体注入影响的裂隙岩石的流体力学行为(参考文献 1)。该模型可以真实地捕捉裂隙多孔介质中的压力扩散,完整岩石中由脆性和疲劳引起的损伤以及裂隙结构对流体力学过程的重要影响(图2)。该模型还使我们能够直观地查看裂隙岩石中损伤、应力和压力场的详细演变,并进一步研究多孔弹性对驱动系统中新损伤传播的基本控制(图 3)。根据模拟结果,我们还可以分析由完整岩石脆性破坏和/或天然裂隙摩擦滑动引起的诱发地震活动的时空演变(图4)。
图2 注液过程中裂隙岩石的压力演化与损伤扩展。
一个 ×ばつ3 的图像网格,显示了破裂岩石局部区域的顶行损伤、中行应力比和底行流体压力
图3 查看(a)损伤的分布情况;(b)应力比(即局部最大主应力与局部最小主应力的比值);(c)裂隙岩石局部区域的流体压力(通过高度表达式显示)。
一个 ×ばつ4 的图像网格,显示了不同时间戳下裂隙岩中诱发地震事件的空间分布和演化
图4 低、高裂隙密度分别为χ = 0.5 和 1.5 的裂隙岩石中诱发地震事件的空间分布和演化规律。
该模型也可用于模拟裂隙岩石中开挖引起的扰动以及由此产生的瞬态流体力学行为(参考文献4)。我们捕获了由于挖掘(时间 t = 0-0.1 小时)和随后的排水(时间 t= 0.1~20 h) 过程导致的显著的压力变化和扩散以及应力变化和损伤演变(图5)。我们通过对 Biot 系数进行敏感性分析来说明流体力学耦合的重要作用。结果表明,当Biot系数越高(或者说耦合越强)时,开挖引起的孔隙弹性压力场越不均匀,岩石损伤和破裂位移也越大。开挖和排水过程都会诱发与岩石基质的脆性损坏和/或天然裂隙的摩擦滑动相关的地震事件(图 6)。
图5 裂隙岩石在开挖过程中及开挖后的压力、应力和损伤演化。
一个 ×ばつ2 的图像网格显示了断裂岩石中地震事件的空间分布,左边的图像显示了挖掘过程中的岩石,右边的图像显示了排水过程中的岩石
图6 在开挖(左图)和排水(右图)阶段,具有不同 Biot 系数 α 的裂隙岩石的地震活动的空间分布。
除了上述的流体力学模型外,我们还开发了一个完全耦合的热流体力学模型来模拟裂隙性地热储层在长期水循环和产热过程中的性能(参考文献3)。
雷庆华博士,讲师,瑞士 ETH Zürich 地球科学系高级科学家。他拥有中国同济大学土木工程学士学位(2009 年)和硕士学位(2012 年),以及英国帝国理工学院岩石力学博士学位(2016 年)。雷博士是国际岩石力学和岩石工程学会(ISRM)的 Rocha 奖章,以及 NGW Cook 博士论文奖和美国岩石力学协会(ARMA)的岩石力学研究奖的获得者。雷博士的研究兴趣包括岩石力学、耦合过程、断裂表征、多相流、地震波、诱发地震活动和边坡稳定性。他是断裂岩石热-水力-机械-化学过程 ISRM 委员会的秘书长,ARMA 未来领导人,以及 ARMA 地下储存和利用技术委员会的创始成员。
Autodesk、Autodesk 徽标、AutoCAD 和 DXF 是 Autodesk, Inc.和/或其附属公司和/或关联公司在美国和/或其他国家的注册商标或商标。
Rhinoceros 是 TLM, Inc. 的注册商标。DBA Robert McNeel & Associates Corporation。
MATLAB 是 The MathWorks, Inc. 的注册商标。
可以详细介绍一下,建模和物理场建模过程么?或者分享一下案例,感谢
本博客暂时还没有相关案例,可以参考官网其他案例:
http://cn.comsol.com/model/discrete-fracture-691
http://cn.comsol.com/model/flow-in-a-fractured-reservoir-90421
雷老师 您好!想学习学习您的建模方法 有基础案例分享吗
雷老师 您好! 想请问您第一个案例中是怎么确保生成的损伤区为线的
您好,几何构建过程如前面的"建模过程"步骤1 所述,需要在几何上将损伤表示为线/折线的形式。参考案例:
http://cn.comsol.com/model/discrete-fracture-691
http://cn.comsol.com/model/flow-in-a-fractured-reservoir-90421
雷老师,您好,请问在固体力学接口您是如何实现对裂隙进行赋值的?感谢!
您好,"多孔弹性"接口可实现固体力学和流体方程之间的直接耦合,您可以在固体力学接口设置您感兴趣的内容。
雷老师您好,想请问您是怎么在相邻单元之间嵌入联合单元的 谢谢!
您好,可以对裂隙几何进行正常的网格划分。
您好 ,您的意思是裂隙在几何阶段就不是边界,是一个宽度远小于长度的矩阵域么
在裂隙边界条件情况下,"联合单元"可以理解为边界条件(边界方程)。 没有创建矩形域网格单元。
相邻单元之间如何嵌入联合单元
您好,您的意思是裂隙几何网格破分生成的自由三角形网格就是已经嵌入联合单元了吗?
老师,请问一下这个模型在后面会分享详细的建模和所用物理场吗?可以通过培训看到这个模型的详细介绍吗?谢谢老师
您好,该博客为用户经验分享,再之后的培训中并不会看到该模型的详细介绍。但COMSOL官方会定期举办岩土相关的研讨会,建议您关注官网相关信息:http://cn.comsol.com/events
官网上的其他相关案例:
http://cn.comsol.com/model/discrete-fracture-691
http://cn.comsol.com/model/flow-in-a-fractured-reservoir-90421
请问流体源项的单位kg/(m3.s),怎么解释呢,和流量m3/s怎么联系呢
可以分享一下案例吗?谢谢
您好,该博客为用户经验分享,在之后的培训中并不会看到该模型的详细介绍。但 COMSOL 官方会定期举办岩土相关的研讨会,建议您关注官网相关信息:http://cn.comsol.com/events
官网上的其他相关案例:
http://cn.comsol.com/model/discrete-fracture-691
http://cn.comsol.com/model/flow-in-a-fractured-reservoir-90421
您好,请问您这个裂纹的路径是提前预制的吗?可以实现裂缝的自由扩展吗?就像扩展有限元那样。
您好,三维那个随机断裂的方案中提到的插件可以分享一下吗?
您好,该插件已在软件中提供,开发工具-插件库。
请问可以在土体中建立上述示例2的模型吗?
您好,示例2中展示了裂隙岩石开挖过程的中压力等属性变化,您所述的土体与裂隙岩石材质类似,可以采用相似的建模方法实现。
您好,我想问下,如果做煤层开挖,在固体力学里加土壤塑性还是岩石还需要加损伤吗,土壤塑性经常找不到解是怎么回事,我用的是剪切模量和体积模量
您好,损伤一般用来仿真裂隙或损伤区域的延申或开展,仿真开挖过程一般不需要考虑损伤,您可以参考案例:http://cn.comsol.com/model/deep-excavation-9747。
您好,请问裂隙的刚度和开度的变化怎么实现呢?
您好,在本篇博客中,裂隙刚度的变化可以根据相应的公式通过定义变量去实现,开度可以通过定义它和裂隙的法向刚度的关系以及考虑其受剪膨胀来进行表征。
请问用固体力学模块的弹性薄层来表征裂隙的力学性质可以达到文献中效果吗
您好,如果是考虑裂缝存在法向刚度以及剪切刚度,可以考虑使用弹性薄层来对裂纹的相关属性进行定义。
请问一下文献是采用哪种方式实现在固体力学部分对裂隙法向、切向位移、开度、刚度等参数设置以及实现裂缝应力集中呢?在固体力学模块没有找到相关接口涉及裂缝开度等参数的定义。
我看COMSOL中的损伤都是拉损伤,请问如何如何同时定义岩土材料的剪切损伤和拉伸损伤呢,譬如剪切损伤阈值为摩尔库伦准则,拉伸损伤阈值为最大拉应力准则
您好,不管是拉伸损伤还是剪切或压缩损伤,都可以通过标量损伤下的等效应力的不同组合方式去考虑。默认的几种组合方式主要考虑三个主应力的影响和组合,如果有一些单独考虑剪应力的组合方式,这里也可以自定义修改。
您好,在等效应变下有一个激活压缩的勾选栏,请问如果勾选了就是变成以压缩损伤作为驱动的损伤了吗?如果想相同获取拉伸与压缩损伤是否需要添加两个损伤节点?
文中提到"我们定义了岩石基质和裂隙的材料特性和本构方程。一些岩石/裂隙特性,例如孔隙率、储水和渗透率,被定义为局部应力/压力状态的函数,以实现间接耦合。"这里岩石基质和裂缝是采用的两个压力变量吗?
您好,文中提到"我们使用 COMSOL Multiphysics 中的固体力学 和达西定律 接口对裂隙介质中的流体力学过程进行了模拟",因此模型仅采用了达西定律一个流动接口来求解流体流动,那么则是在达西定律接口内选择了’裂隙‘特征来模拟裂缝中的流动,这意味着岩石基质和裂缝采用的是一个压力变量,即达西定律的因变量。
请问这种情况下,构建的裂隙可以被认为是离散裂隙介质吗?
请问怎样在comsol中实现像参考案例1那样水力压裂后的裂纹沿着线拓展,固体力学场里需要设置什么呢,后处理上有什么技巧吗?
您好,可以尝试通过固体力学中的"损伤"节点来实现(右键线弹性材料–损伤),并且在后处理中查看"solid.dmg"损伤因子变量。
请问现在有本博客的相关案例吗
您好,本博客暂时还没有相关案例,您可以先参考官网其他案例:
http://cn.comsol.com/model/discrete-fracture-691
http://cn.comsol.com/model/flow-in-a-fractured-reservoir-90421
您好,请问为什么我的损伤是范围伤害,不是线条,而且通过固体力学中的"损伤"节点来实现(右键线弹性材料–损伤),并且在后处理中查看"solid.dmg"损伤因子变量。得到的也是范围性的,不能观察裂缝的扩展路径
您好,根据您的描述这种情况可能是由裂纹正则化(网格)导致。如果采用的是标量损伤,是需要用足够细的网格去合理解析网格(裂缝带)的;如果是采用的相场损伤,相场长度尺寸也应当取到合理值使其能表征裂纹的局部弥散。如果还有其他问题,也可以联系技术支持:http://cn.comsol.com/support
请问示例一是用哪几个模块做的
您好,压裂仿真可能主要会涉及到固体力学以及达西定律接口,二者可以通过预置的多孔弹性接口实现多物理场耦合。
请问固体力学部分的裂缝变形(法向和切向位移)是通过哪个接口设置以实现文献中裂缝尖端的应力集中呢?以及如何在固体力学中实现裂缝开度的设置与变化呢,固体力学中并没有类似达西定律一样可以对裂隙开度进行设置?
您好,如果是考虑裂缝不是完全中空,即存在一定的法向刚度以及剪切刚度,可以考虑使用"薄层"来对裂纹的相关属性进行定义;而如果认为裂纹内部中空,已不具备法向和切向的刚度,则可以尝试使用"裂纹"来进行表征。以上二者特征均可以在固体力学物理场中找到。而关于裂纹的法向位移以及剪切滑移,也都可以在这二者的内置变量中进行查看和调用。
您所说的中空是指裂缝中无流体存在的情况呢,还是指裂缝中没有杂质(如矿物)填充的情况的?您所说的"薄层"是"弹性薄层"还是"薄层"接口呢?
请问如何定义裂隙的法向刚度和剪切刚度呢?薄层节点也没有看到定义这两个属性的
兄弟,你解决了吗?我也没有找到
技术专家您好,关于边界条件"薄层"或"裂隙"或者其他边界条件设置裂隙的变形性质,在"薄层"中的设定"单位长度的力表示为伸长部分的函数"或者"单位长度的弹簧常数"输入接口似乎不能直接定义应力-应变关系,因为法向应变受法向力、切向力共同影响,因此 似乎不易进行转换为"单位长度的力表示为伸长部分的函数",而且也不易变换出应力与应变的比值以导入"单位长度的弹簧常数"。那么如何实现这个本构关系的定义?定义材料属性吗?但"材料数据"中要求给出"弹性模量和泊松比"。另外,我还查阅了"Crack"边界条件,它可以用于这一设置吗?我注意到其中还有一个问题,就是"Crack"中似乎不能设置裂隙初始宽度。请不吝赐教!
请问图3(c),图4,图6是COMSOL出图吗,如何调用呢?
您好,图3(c)可以尝试在表面绘图中添加高度表达式来实现,而另外两个也可以尝试通过二维绘图组中的更多点图中的"点"来实现。
您好,在等效应变下有一个激活压缩的勾选栏,请问如果勾选了就是变成以压缩损伤作为驱动的损伤了吗?如果想相同获取拉伸与压缩损伤是否需要添加两个损伤节点?
您好,对于标量损伤模型来说,如果想同时考虑压缩和拉伸损伤:推荐首先可以在损伤演化曲线中,使用分段函数或if条件语句,针对拉伸和压缩分别定义其损伤(软化)曲线;之后采用同样的方式,在等效应变中分别定义等效拉伸应变和压缩应变即可。
麻烦具体说一下如何体现损伤,此模拟能否体现裂缝渗透率变化
您好,损伤仿真一般通过固体力学中的"损伤"节点来实现(右键线弹性材料–损伤),并且在后处理中查看"solid.dmg"损伤因子变量,通过该标量场的演化来表征损伤的开展。另外,您可以尝试把渗透率和裂缝开度定义为何其有关的变量,就能够实现渗透率的动态演化。
您好,如果我想模拟地下气驱水的过程,需要考虑注气时裂缝的扩展与损伤,剪切滑移,貌似与这篇文章有些相似,可以具体说一下在固体力学和达西定律下都用哪些接口吗
您好,由于6.3版本删除了达西定律物理场中裂隙的源项,如果我想加裂隙源项该如何操作呢,谢谢
可通过裂隙-点-质量源进行添加。
您好,请问这个裂隙的剪切位移怎么表达呢?
您好,具有法向和切向刚度的裂隙,可以尝试使用(弹性)薄层进行建模,之后使用薄层的剪切位移即可表征裂隙的剪切滑移。
技术专家您好,关于边界条件"薄层"或"裂隙"或者其他边界条件设置裂隙的变形性质,在"薄层"中的设定"单位长度的力表示为伸长部分的函数"或者"单位长度的弹簧常数"输入接口似乎不能直接定义应力-应变关系,因为法向应变受法向力、切向力共同影响,因此 似乎不易进行转换为"单位长度的力表示为伸长部分的函数",而且也不易变换出应力与应变的比值以导入"单位长度的弹簧常数"。那么如何实现这个本构关系的定义?定义材料属性吗?但"材料数据"中要求给出"弹性模量和泊松比"。另外,我还查阅了"Crack"边界条件,它可以用于这一设置吗?我注意到其中还有一个问题,就是"Crack"中似乎不能设置裂隙初始宽度。请不吝赐教!
请问,本案例中所用的裂纹损伤模型一般用于准脆性或者混凝土,那么什么情况下可以用于土壤或岩层的裂纹延展?如果不可以应该如何实现呢?
您好,博客中用到的标量损伤模型也可以表征土壤以及岩石的拉伸或压缩破坏,除此之外,如果想考虑(例如摩尔库伦准则驱动的)剪切破坏,也可以通过自定义损伤模型的方式来实现。
技术专家您好,抱歉我才意识到单独留言会更适合。关于边界条件"薄层"或"裂隙"或者其他边界条件设置裂隙的变形性质,在"薄层"中的设定"单位长度的力表示为伸长部分的函数"或者"单位长度的弹簧常数"输入接口似乎不能直接定义应力-应变关系,因为法向应变受法向力、切向力共同影响,因此 似乎不易进行转换为"单位长度的力表示为伸长部分的函数",而且也不易变换出应力与应变的比值以导入"单位长度的弹簧常数"。那么如何实现这个本构关系的定义?定义材料属性吗?但"材料数据"中要求给出"弹性模量和泊松比"。另外,我还查阅了"Crack"边界条件,它可以用于这一设置吗?我注意到其中还有一个问题,就是"Crack"中似乎不能设置裂隙初始宽度。请不吝赐教!
您好,博客中的案例使用的是(弹性)薄层节点来进行设置,其中法向和切向的位移和其对应方向的作用力线性相关。如果您认为"法向应变受法向力、切向力共同影响",可以参考博客中参考文献的方法,将对应的法向刚度通过自定义变量的形式完成定义。
专家您好,我看你在给其他人的回复中说使用薄层接口来模拟雷老师案例中所说的裂缝法向和切向位移行为。我目前使用的是comsol6.2版本,在固体力学下的薄层父接口下有”固体,膜,弹簧”三个选项,还对应有四种材料模型”线弹性材料,超弹性材料,非线性弹性材料,弹簧材料”,请问应该分别选用哪种更好呢?另外,我在这些选项下面都没有找到直接定义法向和切向刚度的地方,弹簧材料下倒是可以定义法向刚度但是无法定义切向。
三种薄层类型中,固体是最全的,适合一般情况,膜假设只有切向应变,而弹簧假设只有法向应变;考虑裂隙的法向刚度和切向刚度可以采用"固体"和"弹簧"。材料类型一般使用线弹性材料即可,如果您只知道刚度值,您还需要结合厚度转换为对应的杨氏模量等值;如果您用弹簧,可以输入一个K矩阵,能定义法向和切向刚度。
您好,之前 Ge 工程师提到 "结合厚度转换为对应的杨氏模量等值",是依据裂缝刚度来设置薄层的厚度以及材料属性(比如杨氏模量、泊松比与密度)吗?关于 "厚度转换为杨氏模量" 的具体转换关系,能麻烦您帮忙解答一下。