如何运用COMSOL Multiphysics来模拟水力压裂和岩石损伤

引言  

为了开采地层中的能源，需要对低渗透性的储层进行压裂，从而形成裂隙网络，便于能源开采。  

目前关于水力压裂的研究很多，主要包括试验和数值模拟，大多数的研究都是将二者结合。今天主要分享一下如何运用有限元软件COMSOL Multiphysics来模拟水力压裂。  

岩石水力压裂（Song et al. 2019, Joule）

在以往的研究中，用COMSOL模拟水力压裂都需要借助MATLAB来传递数据计算损伤。但这样可想而知效率不高，而且麻烦，对一般人来说不易掌握。经过不断的探索，后面我发现其实可以省去MATLAB这个步骤，直接在COMSOL里进行。  

建模过程和心得  

在COMSOL里建模，理解物理过程很重要。水力压裂的主要物理过程就是流体注入，压力增长会对岩石产生损伤，比如说弹性模量。而在COMSOL里直接对弹性模量赋损伤值，则会导致错误：检测到循环变量。其实原因也很简单，因为计算应力应变需要通过弹性模量来计算，而损伤变量中又含有应力应变，所以成了死循环。  

首先看通常所采用的本构模型。最基本的压裂方式是采用水来进行压裂，然而在深部地层中必须考虑热的影响。所以模型多采用水-力-热-损伤（THMD）模型，这样还可以研究不同流体对损伤的影响。在COMSOl里的计算模块主要包含达西定律，固体力学和多孔介质传热。  

COMSOL模块选择  

这里的话只给出三个主要的耦合方程，第一个就是力学方程。在这个方程里一方面需要考虑流体的影响（左边第三项），另一方面需要考虑温度对力学的影响（左边第四项）：  

第二个方程就是和流体压力扩散相关的方程，主要的耦合还是要考虑固体变形对流体的影响（左边第三项）以及温度的影响（左边最后一项）：  

第三个方程就是能量守恒方程，同样需要与压力和温度耦合：  

只有将上面三个方程完全耦合才能算全耦合。  

我看了很多损伤和水力压裂相关的文章，里面的方程可能存在错误，即便是发在Energy这些期刊上（看看我给的方程正不正确）。需要注意的是看公式里耦合项的正负号是否正确，以及正负号的定义，搞错了结果就是错的。  

当损伤产生时，会对孔隙度、渗透率和导热系数等参数产生影响，具体就看你模拟哪些参数了，下面主要列举常见的几个关系：  

孔隙度变化  

渗透率变化  

热传导系数变化  

损伤的本构是采用下面这个方程，这个方程在很多论文里都有采用。一个是压损，一个是拉损。当然还是剪切损伤，具体也可以定义，可能会复杂一点。  

拉损和压损定义  

这里损伤变量D是和应变有关，满足不同条件时发生相应的损伤：  

损伤变量定义  

关于各种变量，直接在COMSOL中用变量定义即可。不同物理场的耦合则需要修改底部方程以及利用COSMOL已有的接口。  

物理场设置  

关于岩石的非均匀参数，例如弹性模量和渗透率等，都是通过MATLAB用Weibull分布生成，然后再导入到COMSOL中。下面给出生成的例子。  

弹性模量分布  

当模型设置好后，最难的一个点就是对参数取值。我看了很多论文，参数基本都不给全。你不给全，那自然是此地无银三百两。反正是模拟嘛，不调参那能叫模拟？习惯了。当然，也可能由于有的参数确实没有测，自己根据别人的文献参考了一下。模拟挺玄乎，损伤也挺玄乎的。  

当模型设置好后，接下来就是漫无止境的调参，这也是COMSOL的精髓。COSMOL这个软件很玄乎，有时候稍微动一下就不收敛。所以对于新手而言，最好先从简单模型入手，这样错了也能知道为什么。  

我目前已经调试了三种不同流体的水力压裂，分别是水、超临界CO2以及液氮。当流体的性质涉及到温压变化时，模型是最难收敛的。所以我水的模型是收敛性最好，滋溜一下就算完了。另外两个模型目前也能计算。  

最近刚好做了一个液氮压裂的模型，模型的尺寸是1*1m，中间一个压裂孔，流体从这里注入。  

压裂模型图  

这个模型要比通常的试验模型大，所以随便计算一下。岩石的初始温度是373K，液氮的注入温度是100K。当液氮注入时会形成一个热冲击，对岩石造成一定损伤。围压分别设置为25和20MPa。所以裂缝会沿垂直于最小主应力的方向发展（看下哪个方向是最小主应力方向）。  

下面直接给出计算结果吧  

损伤分布  

压力分布  

温度分布  

（为了不影响展示结果，未在图片中间加水印，请勿盗图搞销售）

之前算的模型，损伤都没有像液氮这样美的！！！简直太棒了  

当然了这个模型的参数是我随便取的，没有那么严谨。只是尝试一下看能不能计算。  

果然，科研最棒的部分就是出结果那部分！！！