本发明属于能源开采与地下储存，更具体地，涉及一种基于comsol自迭代计算岩石损伤的方法及介质。
背景技术：
1、为了开采非常规能源，如页岩气的开采，这些能源通常存在于低渗透性的储层中，需要通过水力压裂等技术形成裂缝网络。不同作用机制形成的裂缝将影响能源的开采。另一方面，为了减少大气中温室气体的含量，需要将二氧化碳(co2)封存在地层中，大量的co2注入可能导致地层岩石损伤，对盖层的完整性造成影响。流体注入地层的例子还包括废水处置，同样会对岩石造成损伤。通过开展实验模拟岩石损伤具有局限性，因为地层中岩石的损伤过程包含复杂的热-流-固-损伤(t-h-m-d)多场耦合作用。实验过程中难以模拟实际条件，因此需要借助数值计算来模拟岩石损伤。虽然目前有许多模拟岩石损伤的模拟软件和方法，但大都局限于流固耦合，在多场耦合方面没有有限元软件comsol multiphysics那样便于开发。comsol multiphysics是一款强大的多物理场耦合软件，为多场耦合提供了丰富的接口，同时也具有很好的二次开发功能。目前基于comsol模拟岩石损伤需要借助matlab进行应力传递来计算岩石的损伤因子，comsol和matlab的交互使得该方法效率较低，同时不便于多场耦合，以及考虑更复杂的条件。
技术实现思路
1、针对现有技术的以上缺陷或改进需求，本发明提出了一种基于comsol自迭代模拟岩石损伤的方法及介质，以解决现有技术需要借助matlab计算损伤变量的弊端，提高求解效率，便于开展多场耦合研究。
2、为实现上述目的，按照本发明的一个方面，提供了一种基于comsol自迭代计算岩石损伤的方法，包括：
3、确定模拟岩石的力学参数和渗透参数；
4、在matlab中通过weibull分布生成非均匀的材料参数分布，采用插值方式导入comsol；
5、在comsol中定义参数和变量，绘制几何，添加渗流场和力学场，设置边界条件；
6、采用状态变量计算每一步的损伤，对材料的力学参数和渗透参数进行折减；
7、采用delaunay划分方式对模型进行网格划分，同时对注入口处网格局部加密；
8、对模型应力和应变进行初始化，模拟实际注入前的应力状态；
9、添加瞬态计算，选择稳态作为初始值计算模型。
10、在一些可选的实施方案中，岩石力学参数包括弹性模量、泊松比、内摩擦角、抗拉强度和抗压强度，岩石渗透参数包括渗透率和孔隙度。
11、在一些可选的实施方案中，所述在matlab中通过weibull分布生成非均匀的材料参数分布，采用插值方式导入comsol，包括：
12、在matlab中编写weibull分布非均质材料参数的代码，确定weibull分布的非均质指数，以及哪些参数服从weibull分布；
13、将matlab的导出格式设置为excel，第一列为x坐标，第二列为y坐标，第三列为对应坐标的材料参数值；
14、在comsol中选择插值函数插入材料参数值。
15、在一些可选的实施方案中，在comsol中定义的参数包括内摩擦角、模型尺寸和注入压力，变量包括渗透率和孔隙度随压力变化的自定义关系。
16、在一些可选的实施方案中，comsol中渗流场用达西定律进行模拟，力学场用固体力学进行模拟，二者通过多孔弹性模块进行耦合，其中，模型的多孔弹性耦合计算公式为：其中，ρ是流体密度，s是储水系数，p是流体压力，μ是流体动力粘度，k是岩石渗透率，q是源项，α是比奥系数，εvol是体积应变，t是时间，σ’是有效应力，i是单位矩阵，f是体积力。
17、在一些可选的实施方案中，对于矩形岩样，模型左侧和下侧设置为辊支撑，限制垂向变形，上侧和右侧分别设置为相应的应力值。
18、在一些可选的实施方案中，所述采用状态变量计算每一步的损伤，对材料的力学参数和渗透参数进行折减，包括：
19、采用状态变量计算每一个时间步的损伤，并设置为在每一步计算之前计算损伤对材料力学参数的影响，进而对材料的力学参数和渗透参数进行折减，其中，采用的判断岩石损伤的强度计算公式和定义的损伤变量为：
20、
21、其中，f1表示抗拉强度准则，f2表示抗压强度准则，σ1和σ3分别是最大最小主应力，ft是抗拉强度，d为损伤变量，εt和εc分别是岩石的拉应变和压应变，ε1和ε3分别是岩石的最大主应变和最小主应变，φ是内摩擦角，fc是抗压强度，e0是初始弹性模量，e是损伤后的弹性模量。
22、在一些可选的实施方案中，设置稳态初始条件对模型进行初始化，初始条件包括应力和应变，模拟实际注入前的应力状态，不考虑流体注入。
23、在一些可选的实施方案中，所述添加瞬态计算，选择稳态作为初始值计算模型，包括：
24、添加瞬态计算，选择稳态作为初始值，根据材料参数和注入压力大小设置总时间和时间步长，计算模型，并设置计算时显示结果，根据计算结果及时调整模型；如果到达设定的时间未出现损伤或未达到预期目标，调整注入压力或者weibull分布的非均匀指数，重新计算。
25、按照本发明的另一方面，提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一项所述方法的步骤。
26、总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：
27、本发明提供的基于comsol自迭代模拟岩石损伤的方法。该方法无需借助matlab计算损伤变量，计算过程完全在comsol软件中进行，提高了求解效率，能求解任意时间步的损伤，同时也便于开展多场耦合研究，为开展复杂条件下的岩石损伤机理研究提供技术支持。
1.一种基于comsol自迭代计算岩石损伤的方法，其特征在于，包括：
2.根据权利要求1所述的方法，其特征在于，岩石力学参数包括弹性模量、泊松比、内摩擦角、抗拉强度和抗压强度，岩石渗透参数包括渗透率和孔隙度。
3.根据权利要求2所述的方法，其特征在于，所述在matlab中通过weibull分布生成非均匀的材料参数分布，采用插值方式导入comsol，包括：
4.根据权利要求3所述的方法，其特征在于，在comsol中定义的参数包括内摩擦角、模型尺寸和注入压力，变量包括渗透率和孔隙度随压力变化的自定义关系。
5.根据权利要求4所述的方法，其特征在于，comsol中渗流场用达西定律进行模拟，力学场用固体力学进行模拟，二者通过多孔弹性模块进行耦合，其中，模型的多孔弹性耦合计算公式为：其中，ρ是流体密度，s是储水系数，p是流体压力，μ是流体动力粘度，k是岩石渗透率，q是源项，α是比奥系数，εvol是体积应变，t是时间，σ’是有效应力，i是单位矩阵，f是体积力。
6.根据权利要求5所述的方法，其特征在于，对于矩形岩样，模型左侧和下侧设置为辊支撑，限制垂向变形，上侧和右侧分别设置为相应的应力值。
7.根据权利要求6所述的方法，其特征在于，所述采用状态变量计算每一步的损伤，对材料的力学参数和渗透参数进行折减，包括：
8.根据权利要求7所述的方法，其特征在于，设置稳态初始条件对模型进行初始化，初始条件包括应力和应变，模拟实际注入前的应力状态，不考虑流体注入。
9.根据权利要求8所述的方法，其特征在于，所述添加瞬态计算，选择稳态作为初始值计算模型，包括：
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至9任一项所述方法的步骤。
