CN114444230A - 一种超临界co2作用下准脆性材料变形-碎裂的模拟方法 - Google Patents

一种超临界co2作用下准脆性材料变形-碎裂的模拟方法 Download PDF

Info

Publication number
CN114444230A
CN114444230A CN202210210428.3A CN202210210428A CN114444230A CN 114444230 A CN114444230 A CN 114444230A CN 202210210428 A CN202210210428 A CN 202210210428A CN 114444230 A CN114444230 A CN 114444230A
Authority
CN
China
Prior art keywords
fracture
crack
unit
stress
node
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210210428.3A
Other languages
English (en)
Inventor
李�浩
韦婕
姚宏波
贺伟
周忠斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Taiyuan University of Technology
Original Assignee
Taiyuan University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Taiyuan University of Technology filed Critical Taiyuan University of Technology
Priority to CN202210210428.3A priority Critical patent/CN114444230A/zh
Publication of CN114444230A publication Critical patent/CN114444230A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C60/00—Computational materials science, i.e. ICT specially adapted for investigating the physical or chemical properties of materials or phenomena associated with their design, synthesis, processing, characterisation or utilisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computer Hardware Design (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明公开了一种超临界CO2作用下准脆性材料变形‑碎裂的模拟方法，包括以下步骤：统计煤岩体中天然裂隙分布，确定几何参数，编写裂隙生成程序；建立数值计算模型，划分实体单元，在实现天然裂隙及实体单元交界处嵌入裂缝单元并自动生成集合；依次构建反映煤岩基质塑性变形、非贯通裂缝混合型断裂、煤岩完全断裂后块体间接触面分离/挤压/压剪摩擦的力学本构关系；确定煤岩力学参数随超临界CO2作用时间的演化规律，输入材料参数并设置初始条件；设置增量步大小、输出参量、裂缝单元失效及接触对激活时刻的状态变量。本发明为CO2地质封存条件下，盖层受CO2腐蚀、复杂应力耦合作用下其密封性演化提供了有效的模拟工具。

Description

一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法
技术领域
本发明涉及CO2地质封存、碳减排技术领域，具体涉及一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法。
背景技术
CCUS(Carbon Capture,Utilization and Storage，碳捕集、利用与封存)技术是减少化石能源发电和工业生产过程中CO2排放的关键技术，也是我国实现碳中和的兜底技术。CO2地质封存是CCUS技术的核心组成部分，决定了CCUS技术的发展潜力和发展方向。CO2长期安全封存的关键是在限定地质储层环境中不发生泄漏，核心要素是CO2长期作用与深部应力耦合作用下盖层岩石的变形-碎裂力学性能变化特征。在深部煤层中，CO2易达到超临界状态条件(7.38MPa，31.4℃)，这种腐蚀性流体会对盖层岩石性状及其力学性能产生重大影响，采用合理的数值模拟方法预测超临界CO2-应力耦合作用下盖层变形-碎裂发育高度，对于定量分析CO2泄漏规律、优化CO2注入参数、封存位置等具有重大意义。
近年来，国内外学者主要采用有限元和离散元方法开展CO2地质封存、煤岩等准脆性介质的变形与破裂模拟研究。基于有限元方法，采用弹塑性力学理论、损伤力学理论，研究CO2注入地层引起盖层的应力、位移、塑性区、损伤区等的影响范围。基于离散元方法，将煤岩简化为砌体、球体等基本部件组成的集合体，完整煤岩采用刚体或弹性体，裂缝部分多采用库伦滑移模型，进而数值计算复杂应力状态下煤岩的碎裂过程。但是，盖层中裂缝是随机分布的，在地应力、CO2注入压力联合影响下，上述完整煤岩体、裂缝往往受到复杂应力作用，导致其出现塑性变形、混合型断裂模式、分离/挤压/剪切摩擦等复杂响应过程；更重要的是，超临界CO2会腐蚀岩石内的石英、方解石、黏土矿物等成分，导致其强度、断裂能、剪切力学性能降低。忽视这二个因素将会严重影响CO2地质封存数值模拟结果的准确性。
发明内容
为解决上述问题，本发明提供了一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法，可分析深部应力-超临界CO2长期腐蚀耦合作用下盖层变形-碎裂机理，CO2沿盖层及其内裂缝泄漏规律，为CO2地质封存参数优化提供有效模拟手段。
为实现上述目的，本发明采取的技术方案为：
一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法，包括以下步骤：
S1、统计煤岩体中天然裂隙分布，确定几何参数，编写裂隙生成程序；
S2、建立数值计算模型，划分实体单元，在实现天然裂隙及实体单元交界处嵌入裂缝单元并自动生成集合；
S3、依次构建反映煤岩基质塑性变形、非贯通裂缝混合型断裂、煤岩完全断裂后块体间接触面分离/挤压/压剪摩擦的力学本构关系；编制有限元-离散元FDEM计算程序，实现煤岩塑性变形-混合型断裂-碎裂摩擦全过程；
S4、实验确定煤岩力学参数随超临界CO2作用时间的演化规律，输入材料参数并设置初始条件；
S5、设置增量步大小、输出参量、裂缝单元失效及接触对激活时刻的状态变量，分别通过有限元、离散元计算单元失效前后模型响应；
S6、分析超临界CO2作用不同时间下，受载煤岩试件变形-碎裂结果。
进一步地，所述步骤S1中，采用CT扫描方法确定煤岩中主要裂缝的分布，统计不同方位角裂隙出现的频数，统计相邻裂缝的平均间距及其方差；在此基础上编制主要裂缝自动生成程序，具体为：根据扫描确定的主要裂缝平均间距及其方差，在MATLAB中随机生成离散点，确定离散点的坐标，将散点连接成Delaunay三角形，由三角形相邻边的中垂线围成Voronoi多面体，生成多组Voronoi多面体，再优选出与CT扫描所得主要裂缝方位角结果接近的Voronoi多面体集合；由上述Voronoi多面体的各个面模拟煤岩中主要裂缝空间分布。
进一步地，所述步骤S2包括如下步骤：
S2.1、建立数值计算模型，划分实体单元网格；
在步骤S1的基础上，进一步根据数值模拟目的确定研究对象形状及尺寸，采用“切割”、“倒圆角”等工具将Voronoi多面体模型截割成所需的数值模型；
网格划分前，依据下式确定单元最大尺寸L：
L＝9πEFe/[32(1-μ2)σp]
式中，E为弹性模量；μ为泊松比；Fe为断裂能；σp为抗拉强度；
通过Mesh模块“Seed edge”对数值模型布种，种子间距为L，进而通过“Seed part”命令对数值模型划分实体单元网格；
S2.2、在天然裂隙及实体单元交界处嵌入裂缝单元，并自动划分集合。
(1)更新实体单元节点编号；
网格划分后，通过CAE生成“.inp文件”，将文件内关键字“*nset”和“*elset”下方、“*asssembly”上方所有数字复制在Excel中，找到最大节点编号nmax、最大单元编号emax；
保持所有节点坐标不变，将第ni(0<i≤nmax)个节点复制a次，同时保持单元编号最小的节点编号ni,min不变，将其余单元的节点编号增加至最大节点编号的[(a-1)×10+ni,min]，其中，a为某个节点共用单元的个数；
(2)创建裂缝单元节点编号、单元编号；
若以四边形实体单元为例，记上述重新编号后的节点编号为N01,N02,N03,N04；共用第n0i个节点的相邻单元的编号为e01,e02,e03,e04。将新的节点编号和对应单元组合成有序数组：
Array＝(e01,e02,e03,e04,N01,N02,N03,N04)
将N01～N04逆序排列，即可得到四边形实体单元e01,e02,e03,e04之间裂缝单元的节点编号；
而裂缝单元最小的编号由实体单元最大编号所在数位+1得到的10的整数幂次方确定；例如实体单元最大节点编号为756，裂缝单元最小编号取值为1000。
(3)在主要裂隙处嵌入裂缝单元；
在“.inp文件”中关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohelemAll”；按照第一列为裂缝单元编号、2～5列为节点号的顺序复制进相应的ei和Ni数字，从而在所有主要裂隙、所有实体单元边界处嵌入零厚度的裂缝单元；
(4)自动划分裂缝单元集合；
复制裂缝单元编号ei及其节点坐标Ni(xi,yi,zi)至Excel中，采用内积法求解所有裂缝单元的方向向量dcoh＝(d1,d2,d3)，以及Voronoi多面体各个面的方向向量dip，若二者共线，再通过空间向量方法计算裂缝单元上任意点到Voronoi面的距离，若为0，则判定裂缝单元位于Voronoi多面体的面上；将所有符合上述三个条件的裂缝单元编号、节点编号复制在.inp文件中，在关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohElem01”，从而将所有主要裂隙上的裂缝单元组成一个名称为“cohelem01”的集合；
在此基础上，通过布尔运算获得所有实体单元边界处的裂缝单元集合“cohelem02”。
进一步地，所述步骤S3包括如下步骤：
S3.1、完整煤岩体弹塑性变形；
弹塑性变形阶段，各要素服从广义胡克定律：
σ＝E0:(ε-εp)
其中，σ为应力，E0为弹性刚度矩阵，ε为总应变。塑性应变
Figure BDA0003533011450000058
由屈服函数、塑性势函数确定，屈服函数F表达式为：
Figure BDA0003533011450000051
其中，F为屈服函数，<g>是Macaulay符号，
Figure BDA0003533011450000052
Figure BDA0003533011450000053
p＝trace(σ)/3，γ＝3(1-Kc)/(2Kc-1)，q＝[2(S:S)/3]1/2,S＝σ+pI式中，σi,max(i＝1,2,3)为最大有效主应力，
Figure BDA0003533011450000054
为双轴、单轴压缩屈服应力之比，σt与σc(εp)分别为拉压应力；Kc为控制屈服面在偏平面形状的参数，岩石类材料取值为0.667；
G＝[(δσttanψ)2+q2]-p tanψ
其中，δ为q-p平面内塑性势函数子午线在张拉区段的尖点曲率，一般取值为0.1；ψ为剪胀角；
S3.2、非贯通裂缝混合型断裂；
岩石断裂过程分为弹性变形和混合型断裂两个阶段；
弹性阶段本构关系为：
σc＝D0,cεc
一旦达到如下条件，进入韧性断裂过程：
Figure BDA0003533011450000055
其中，D0,c是弹性刚度矩阵，σc,n,σc,s,σc,t
Figure BDA0003533011450000056
是法向和两个切向方向的应力(峰值应力)，εc是应变；
为推导混合型韧性断裂本构方程，首先建立拉、剪、以及混合模式下断裂能Gc,n,Gc,S(Gc,S＝Gc,s+Gc,t，Gc,s与Gc,t为二个切向方向的断裂能),Gc,m表达式为：
Figure BDA0003533011450000057
其中，j表示n,S,m(即张拉、剪切、混合型)断裂模式；
断裂能混合比表达式为：
ξ＝Gc,S/(Gc,n+Gc,S)
煤岩完全断裂的判据为：
Figure BDA0003533011450000061
其中，
Figure BDA0003533011450000062
χ分别为完全断裂时的张拉、剪切型断裂能及材料参数；
拉、剪、混合型断裂模式下，各要素本构方程为：
Figure BDA0003533011450000063
其中，dc,j，Dc,j，σc,j，Sc,j分别为损伤变量、弹性模量、应力、塑性位移和总位移，且有：
Figure BDA0003533011450000064
其中，Sc,n与Sc,S分别为纯拉、剪条件下的位移；
假设混合型断裂模式下的塑性位移与剪切塑性位移相等，即有
Figure BDA0003533011450000065
其中，未知参数
Figure BDA0003533011450000066
通过以下步骤确定：①采用三点弯曲试验、加卸载直剪试验，获得拉、剪及混合模式下的σc,m-Sc,m曲线、σc,n-Sc,n曲线、σc,S-Sc,S曲线，计算相应的损伤演化曲线及断裂能
Figure BDA0003533011450000067
②通过“煤岩完全断裂判据表达式”对实验所得多组
Figure BDA0003533011450000068
数据进行拟合，得到相同材料参数χ下的断裂能混合比ξ；③在此基础上联立“峰值应力方程”，得到某个ξ下的峰值荷载
Figure BDA0003533011450000069
对应的σc,n与σc,S值，根据弹性表达式计算峰值位移分量Sc,n,Sc,S；④联立式“断裂能表达式”、“本构方程”，得到任意σc,m对应的σc,n与σc,S分量以及相对应的位移分量Sc,n,Sc,S；⑤将混合模式下的应力、位移替换为纯拉、纯剪的力学参数，并由“本构方程”计算σc,n，σc,S和Sc,n,Sc,S条件下的塑性应变、拉/剪损伤变量；⑦由“本构方程”计算得到得到任意混合比下的dc,m-Sc,m关系；由此确定混合型韧性断裂本构方程；
S3.3、贯通裂缝或块体的分离/挤压/剪切摩擦；
一旦断裂能满足断裂能表达式时，岩石完全断裂，在数值计算中删除该断裂面上的零厚度裂缝单元，消除其力学影响；
在此基础上，首先建立块体间分离、压缩、剪切摩擦关系判据：
①分离：相邻实体单元的任意二个节点间距l＞0，发生分离；
②挤压：l＝0且二者间存在压应力；
③剪切摩擦：若l＝0且二者间存在沿结构面的剪应力；
若相邻块体分离(l＞0)，二者间不存在相互作用力，块体运动服从牛顿第二定律；
若发生挤压，则结构面法向方向的力学本构关系为：
σn＝DnNmaxn/(Nmax-n)
其中，σn为压缩应力；Nmax为结构面最大闭合量，根据三维形貌扫描确定；n为结构面闭合量，变量；Dn为结构面法向模量，取值为其相邻岩块参数；
若发生剪切摩擦，粗糙结构面压剪摩擦本构方程：
Figure BDA0003533011450000071
其中，σS(σS,p)为(峰值)剪应力；DS为剪切刚度；S(sp)为(峰值)剪切位移；参数p由σS,p、Sp、残余剪应力σS,r和残余剪位移Sr确定；
步骤S3.4中，有限元FEM和离散元DEM耦合求解；
以被删除的零厚度裂缝单元为边界，将其内部分视作一个整体，根据部步骤S3.3方法，通过DEM方法求解所有离散块体的节点力；将此节点力作为边界条件，根据部步骤S3.1方法，通过FEM方法求解其内部完整块体的节点力；
通过实体单元与裂缝单元的共享节点传递力场数据，在此基础上，根据步骤S3.1计算实体单元的塑性应变、应力σ响应；根据步骤S3.2计算裂缝单元的损伤、位移、应力σc响应；若有断裂能
Figure BDA0003533011450000072
同时|σ-σc|≤10-5，则直接将实体单元、裂缝单元应力、位移等响应传递至下一计算步骤；
若断裂能
Figure BDA0003533011450000081
表明裂缝单元完全断裂，此时删除零厚度裂缝单元，同时根据步骤S3.3计算煤岩块体间的分离/挤压/剪切摩擦；
若同时|σ-σc|＞10-5，则将附加节点力取值为|σ-σc|，根据步骤S3.1和S3.2重新计算实体单元、裂缝单元的节点力，直至|σ-σc|≤10-5成立；而后将实体单元、裂缝单元应力、位移等响应传递至下一计算步骤，重复上述计算过程。
进一步地，所述步骤S4中，将煤岩试件在超临界二氧化碳(ScCO2)制取装置中浸泡0～60d，而后针对完整煤岩试件、含V型缺口的圆柱形试件、含粗糙贯通裂缝试件，分别采用“三轴压力试验机”、“电子万能试验机”、“多功能岩石三轴动态剪切试验机”，获取不同ScCO2浸泡时间下各试件的力学参数；
基于以上实验结果，在“Property”模块下设置1个场变量，场变量数值variables＝试件在ScCO2中的浸泡天数，同时设置相对应的“弹性模量”、“泊松比”、“拉/剪峰值应力”、“摩擦位移”等参数；根据实验条件或工程条件，在“Load”模块下设置相应的应力边界、位移边界、初始应力场等。
进一步地，所述步骤S5中，在“Mesh”模块的“verify mesh”中，检查实体单元的最小稳定增量步，由此确定离散元DEM的最小增量时间步；
在“Mesh”模块下的Element type下设置裂缝单元失效及接触对激活时的变量；
在“Step”模块下的create field output中设置需要输出的变量，包括应力、应变、塑性应变、损伤变量、断裂能、摩擦应力、摩擦位移等；
在此基础上，提交数值计算。
本发明具有以下有益效果：
1)本发明为CO2地质封存条件下，盖层受CO2腐蚀、复杂应力耦合作用下其密封性演化提供了有效的模拟工具。采用本发明所提供的力学本构关系，以及ScCO2浸泡后煤岩力学参数获取办法及参数输入方法，可以实现任意ScCO2腐蚀时间下，煤岩等准脆性材料在荷载条件下的弹塑性变形→损伤→断裂→块体间分离/挤压/剪切摩擦全过程数值模拟。
2)该方法根据断裂能实现裂缝单元删除，在此基础上划定显示时间积分的离散元以及隐式时间积分的有限元计算的空间区域，可以大大提高数值解的收敛性，同时加快计算效率。
3)全局嵌入裂缝单元方法，可实现天然大尺度裂缝、煤岩块体内部微裂缝自动分组，结合有限离散元FDEM方法，最终实现CO2腐蚀、复杂应力耦合作用下裂缝在煤岩中的任意扩展、随机碎裂。
4)本发明可用于分析不同CO2注入压力、不同地质参数影响下CO2储层及盖层的裂缝萌生、生长、扩展演化规律，从而优化CO2地质封存位置、注入压力，实现CO2长期安全封存。助力我国“2030碳达峰、2060碳中和”目标实现。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1为数值模拟方法流程示意图。
图2为CT扫描煤岩实验结果。
图3为数值计算模型及其内贯通、非贯通裂缝分布；
图中；1-主要裂缝处的裂缝单元；2-实体单元间的裂缝单元。
图4为混合型断裂模式的力学响应示意图。
图5为三点弯曲、贯穿剪切断裂力学实验图。
图中：(a)I型断裂试件及其尺寸；(b)II型断裂试件及其尺寸；(c)张拉损伤演化曲线；(d)剪切损伤演化曲线。
图6为粗糙断面三维形貌及其直剪实验图。
图7为为超临界CO2浸泡的数值计算结果；
图中：(a)为超临界CO2浸泡0d的数值计算结果；(b)为超临界CO2浸泡15d的数值计算结果；(c)为超临界CO2浸泡60d的数值计算结果。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。
如图1所示，本发明实施例的一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法，包括以下步骤：
S1、统计煤岩体中天然裂隙分布，确定几何参数，编写裂隙生成程序；具体的，采用CT扫描方法确定煤岩中主要裂缝的分布，统计不同方位角裂隙出现的频数，统计相邻裂缝的平均间距及其方差；在此基础上编制主要裂缝自动生成程序，具体为：根据扫描确定的主要裂缝平均间距及其方差，在MATLAB中随机生成离散点，确定离散点的坐标，将散点连接成Delaunay三角形，由三角形相邻边的中垂线围成Voronoi多面体，生成多组Voronoi多面体，再优选出与CT扫描所得主要裂缝方位角结果接近的Voronoi多面体集合；由上述Voronoi多面体的各个面模拟煤岩中主要裂缝空间分布；
S2、建立数值计算模型，划分实体单元，在实现天然裂隙及实体单元交界处嵌入裂缝单元并自动生成集合；
S2.1、建立数值计算模型，划分实体单元网格；
在步骤S1的基础上，进一步根据数值模拟目的确定研究对象形状及尺寸，采用“切割”、“倒圆角”等工具将Voronoi多面体模型截割成所需的数值模型；
网格划分前，依据下式确定单元最大尺寸L：
L＝9πEFe/[32(1-μ2)σp]
式中，E为弹性模量；μ为泊松比；Fe为断裂能；σp为抗拉强度；
通过Mesh模块“Seed edge”对数值模型布种，种子间距为L，进而通过“Seed part”命令对数值模型划分实体单元网格；
S2.2、在天然裂隙及实体单元交界处嵌入裂缝单元，并自动划分集合。
(1)更新实体单元节点编号；
网格划分后，通过CAE生成“.inp文件”，将文件内关键字“*nset”和“*elset”下方、“*asssembly”上方所有数字复制在Excel中，找到最大节点编号nmax、最大单元编号emax；
保持所有节点坐标不变，将第ni(0<i≤nmax)个节点复制a次，同时保持单元编号最小的节点编号ni,min不变，将其余单元的节点编号增加至最大节点编号的[(a-1)×10+ni,min]，其中，a为某个节点共用单元的个数；
(2)创建裂缝单元节点编号、单元编号；
若以四边形实体单元为例，记上述重新编号后的节点编号为N01,N02,N03,N04；共用第n0i个节点的相邻单元的编号为e01,e02,e03,e04。将新的节点编号和对应单元组合成有序数组：
Array＝(e01,e02,e03,e04,N01,N02,N03,N04)
将N01～N04逆序排列，即可得到四边形实体单元e01,e02,e03,e04之间裂缝单元的节点编号；
而裂缝单元最小的编号由实体单元最大编号所在数位+1得到的10的整数幂次方确定；例如实体单元最大节点编号为756，裂缝单元最小编号取值为1000。
(3)在主要裂隙处嵌入裂缝单元；
在“.inp文件”中关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohelemAll”；按照第一列为裂缝单元编号、2～5列为节点号的顺序复制进相应的ei和Ni数字，从而在所有主要裂隙、所有实体单元边界处嵌入零厚度的裂缝单元；
(4)自动划分裂缝单元集合；
复制裂缝单元编号ei及其节点坐标Ni(xi,yi,zi)至Excel中，采用内积法求解所有裂缝单元的方向向量dcoh＝(d1,d2,d3)，以及Voronoi多面体各个面的方向向量dip，若二者共线，再通过空间向量方法计算裂缝单元上任意点到Voronoi面的距离，若为0，则判定裂缝单元位于Voronoi多面体的面上；将所有符合上述三个条件的裂缝单元编号、节点编号复制在.inp文件中，在关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohElem01”，从而将所有主要裂隙上的裂缝单元组成一个名称为“cohelem01”的集合；
在此基础上，通过布尔运算获得所有实体单元边界处的裂缝单元集合“cohelem02”。
S3、依次构建反映煤岩基质塑性变形、非贯通裂缝混合型断裂、煤岩完全断裂后块体间接触面分离/挤压/压剪摩擦的力学本构关系；编制有限元-离散元FDEM计算程序，实现煤岩塑性变形-混合型断裂-碎裂摩擦全过程；
S3.1、完整煤岩体弹塑性变形；
弹塑性变形阶段，各要素服从广义胡克定律：
σ＝E0:(ε-εp)
其中，σ为应力，E0为弹性刚度矩阵，ε为总应变。塑性应变
Figure BDA0003533011450000125
由屈服函数、塑性势函数确定，屈服函数F表达式为：
Figure BDA0003533011450000121
其中，F为屈服函数，<g>是Macaulay符号，
Figure BDA0003533011450000122
Figure BDA0003533011450000123
p＝trace(σ)/3，γ＝3(1-Kc)/(2Kc-1)，q＝[2(S:S)/3]1/2,S＝σ+pI式中，σi,max(i＝1,2,3)为最大有效主应力，
Figure BDA0003533011450000124
为双轴、单轴压缩屈服应力之比，σt与σc(εp)分别为拉压应力；Kc为控制屈服面在偏平面形状的参数，岩石类材料取值为0.667；
G＝[(δσttanψ)2+q2]-p tanψ
其中，δ为q-p平面内塑性势函数子午线在张拉区段的尖点曲率，一般取值为0.1；ψ为剪胀角；
S3.2、非贯通裂缝混合型断裂；
岩石断裂过程分为弹性变形和混合型断裂两个阶段；
弹性阶段本构关系为：
σc＝D0,cεc
一旦达到如下条件，进入韧性断裂过程：
Figure BDA0003533011450000131
其中，D0,c是弹性刚度矩阵，σc,n,σc,s,σc,t
Figure BDA0003533011450000132
是法向和两个切向方向的应力(峰值应力)，εc是应变；
为推导混合型韧性断裂本构方程，首先建立拉、剪、以及混合模式下断裂能Gc,n,Gc,S(Gc,S＝Gc,s+Gc,t，Gc,s与Gc,t为二个切向方向的断裂能),Gc,m表达式为：
Figure BDA0003533011450000133
其中，j表示n,S,m(即张拉、剪切、混合型)断裂模式；
断裂能混合比表达式为：
ξ＝Gc,S/(Gc,n+Gc,S)
煤岩完全断裂的判据为：
Figure BDA0003533011450000134
其中，
Figure BDA0003533011450000135
χ分别为完全断裂时的张拉、剪切型断裂能及材料参数；
拉、剪、混合型断裂模式下，各要素本构方程为：
Figure BDA0003533011450000136
其中，dc,j，Dc,j，σc,j，Sc,j分别为损伤变量、弹性模量、应力、塑性位移和总位移，且有：
Figure BDA0003533011450000137
其中，Sc,n与Sc,S分别为纯拉、剪条件下的位移；
假设混合型断裂模式下的塑性位移与剪切塑性位移相等，即有
Figure BDA0003533011450000141
其中，未知参数
Figure BDA0003533011450000142
通过以下步骤确定：①采用三点弯曲试验、加卸载直剪试验，获得拉、剪及混合模式下的σc,m-Sc,m曲线、σc,n-Sc,n曲线、σc,S-Sc,S曲线，计算相应的损伤演化曲线及断裂能
Figure BDA0003533011450000143
②通过“煤岩完全断裂判据表达式”对实验所得多组
Figure BDA0003533011450000144
数据进行拟合，得到相同材料参数χ下的断裂能混合比ξ；③在此基础上联立“峰值应力方程”，得到某个ξ下的峰值荷载
Figure BDA0003533011450000145
对应的σc,n与σc,S值，根据弹性表达式计算峰值位移分量Sc,n,Sc,S；④联立式“断裂能表达式”、“本构方程”，得到任意σc,m对应的σc,n与σc,S分量以及相对应的位移分量Sc,n,Sc,S；⑤将混合模式下的应力、位移替换为纯拉、纯剪的力学参数，并由“本构方程”计算σc,n，σc,S和Sc,n,Sc,S条件下的塑性应变、拉/剪损伤变量；⑦由“本构方程”计算得到得到任意混合比下的dc,m-Sc,m关系；由此确定混合型韧性断裂本构方程；
S3.3、贯通裂缝或块体的分离/挤压/剪切摩擦；
一旦断裂能满足断裂能表达式时，岩石完全断裂，在数值计算中删除该断裂面上的零厚度裂缝单元，消除其力学影响；
在此基础上，首先建立块体间分离、压缩、剪切摩擦关系判据：
①分离：相邻实体单元的任意二个节点间距l＞0，发生分离；
②挤压：l＝0且二者间存在压应力；
③剪切摩擦：若l＝0且二者间存在沿结构面的剪应力；
若相邻块体分离(l＞0)，二者间不存在相互作用力，块体运动服从牛顿第二定律；
若发生挤压，则结构面法向方向的力学本构关系为：
σn＝DnNmaxn/(Nmax-n)
其中，σn为压缩应力；Nmax为结构面最大闭合量，根据三维形貌扫描确定；n为结构面闭合量，变量；Dn为结构面法向模量，取值为其相邻岩块参数；
若发生剪切摩擦，粗糙结构面压剪摩擦本构方程：
Figure BDA0003533011450000151
其中，σS(σS,p)为(峰值)剪应力；DS为剪切刚度；S(sp)为(峰值)剪切位移；参数p由σS,p、Sp、残余剪应力σS,r和残余剪位移Sr确定；
步骤S3.4中，有限元FEM和离散元DEM耦合求解；
以被删除的零厚度裂缝单元为边界，将其内部分视作一个整体，根据部步骤S3.3方法，通过DEM方法求解所有离散块体的节点力；将此节点力作为边界条件，根据部步骤S3.1方法，通过FEM方法求解其内部完整块体的节点力；
通过实体单元与裂缝单元的共享节点传递力场数据，在此基础上，根据步骤S3.1计算实体单元的塑性应变、应力σ响应；根据步骤S3.2计算裂缝单元的损伤、位移、应力σc响应；若有断裂能
Figure BDA0003533011450000152
同时|σ-σc|≤10-5，则直接将实体单元、裂缝单元应力、位移等响应传递至下一计算步骤；
若断裂能
Figure BDA0003533011450000153
表明裂缝单元完全断裂，此时删除零厚度裂缝单元，同时根据步骤S3.3计算煤岩块体间的分离/挤压/剪切摩擦；
若同时|σ-σc|＞10-5，则将附加节点力取值为|σ-σc|，根据步骤S3.1和S3.2重新计算实体单元、裂缝单元的节点力，直至|σ-σc|≤10-5成立；而后将实体单元、裂缝单元应力、位移等响应传递至下一计算步骤，重复上述计算过程
S4、实验确定煤岩力学参数随超临界CO2作用时间的演化规律，输入材料参数并设置初始条件；具体的，将煤岩试件在超临界二氧化碳(ScCO2)制取装置中浸泡0～60d，而后针对完整煤岩试件、含V型缺口的圆柱形试件、含粗糙贯通裂缝试件，分别采用“三轴压力试验机”、“电子万能试验机”、“多功能岩石三轴动态剪切试验机”，获取不同ScCO2浸泡时间下各试件的力学参数；基于以上实验结果，在“Property”模块下设置1个场变量，场变量数值variables＝试件在ScCO2中的浸泡天数，同时设置相对应的“弹性模量”、“泊松比”、“拉/剪峰值应力”、“摩擦位移”等参数；根据实验条件或工程条件，在“Load”模块下设置相应的应力边界、位移边界、初始应力场等；
S5、设置增量步大小、输出参量、裂缝单元失效及接触对激活时刻的状态变量，分别通过有限元、离散元计算单元失效前后模型响应；在“Mesh”模块的“verify mesh”中，检查实体单元的最小稳定增量步，由此确定离散元DEM的最小增量时间步；
在“Mesh”模块下的Element type下设置裂缝单元失效及接触对激活时的变量；
在“Step”模块下的create field output中设置需要输出的变量，包括应力、应变、塑性应变、损伤变量、断裂能、摩擦应力、摩擦位移等；
在此基础上，提交数值计算。
S6、分析超临界CO2作用不同时间下，受载煤岩试件变形-碎裂结果。
实施例1
S1.将晋城矿区煤岩试样加工成100mm×100mm×50mm的试件，通过CT扫描获取其内天然裂缝分布，如图2所示。统计天然裂缝间距、与短边方向的夹角等几何参数，编制天然裂缝自动生成程序。
S2.在“Part模块”中，根据S1步骤生成多个Voronoi多面体，以CT扫描结果为依据，选出最接近的Voronoi多面体，并通过切削方法建立直径50mm，高度100mm的圆柱体数值计算模型。计算得到裂缝单元的最大尺寸为6mm，在“Mesh模块”下对数值模型整体划分四面体实体单元网格。而后通过全局嵌入裂缝单元的方法，分别在数值模型的天然裂缝处、实体单元边界处嵌入三维6节点的裂缝单元，并分别组成“cohelem01”和“cohelem02”的裂缝单元集合，如图3所示。
S3.根据煤岩基质塑性变形、非贯通裂缝混合型断裂、煤岩完全断裂后块体间接触面分离/挤压/压剪摩擦的力学本构关系，编制相应的有限元-离散元FDEM的Fortran计算程序，实现煤岩塑性变形-混合型断裂-碎裂摩擦全过程。
S4.将圆柱形试件、含V型缺口的圆柱体试件、含圆环形缺口的圆柱体试件、含粗糙断裂面的立方体试件置于含有超临界CO2的高压反应釜中，静置0d～60d。取出后分别开展单/三轴压缩实验、三点弯曲实验、贯穿剪切实验、剪切摩擦试验、粗糙断面三维形貌扫描，分别获取煤岩块体的弹性模量、泊松比、总应变、塑性应变、三向峰值应力与单轴峰值应力之比、剪胀角，将上述参数赋予实体单元；煤岩非贯通结构面的拉/剪/混合型模量、拉/剪峰值应力、材料参数，将上述参数赋予裂缝单元；粗糙断面天然间距；压剪摩擦过程中的剪切模量、峰值荷载、峰值位移、残余位移等参数，将上述参数赋予裂缝单元失效后激活的接触对。不同超临界CO2浸泡时间对应的参数如表1所示。
表1力学参数
Figure BDA0003533011450000171
S5.在Part模块下生成直径为60mm的圆盘形离散刚体部件，模拟实验过程中压头荷载，并在Assembly模块下复制2个离散刚体部件，设置该部件圆心与圆柱煤岩试件轴线重合，且分别与圆柱上下底面贴合，模拟“压头-试件-底座”系统。在Load模块中设置上部刚体加载速度为0.002mm/min，下部刚体约束全部位移。
根据Mesh模块下verify mesh确定四面体实体单元的最小稳定增量步为10-8，在Step模块中设置一个分析步，其时间增量步为10-8。设置输出变量为应力、应变、塑性应变、断裂能、摩擦位移、摩擦应力、状态变量STATUS等。
S6.在Job模块下新建一个任务，在User subroutine file选项卡中输入编制好的FDEM子程序；在Edit Job选项卡中设置并行计算核心数量为90核，提交计算。
由此得到超临界CO2作用0～60d的数值计算结果，如图7所示。结果表明，超临界CO2作用时间越长，荷载下煤岩碎裂程度越高，盖层长期密封性越低。
本发明与现有受载煤岩破坏模型相比，一方面考虑了超临界CO2对煤岩的腐蚀、强度降低效应；另一方面，考虑了完整煤岩体弹塑性变形、非贯通裂缝混合型断裂、贯通裂缝的分离/挤压/剪切摩擦效应。在此基础上，通过FDEM数值方法模拟了煤岩体变形-碎裂全过程，为超临界CO2长期安全稳定地质封存提供有力方法。
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims (6)

1.一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法，其特征在于：包括以下步骤：
S1、统计煤岩体中天然裂隙分布，确定几何参数，编写裂隙生成程序；
S2、建立数值计算模型，划分实体单元，在实现天然裂隙及实体单元交界处嵌入裂缝单元并自动生成集合；
S3、依次构建反映煤岩基质塑性变形、非贯通裂缝混合型断裂、煤岩完全断裂后块体间接触面分离/挤压/压剪摩擦的力学本构关系；编制有限元-离散元FDEM计算程序，实现煤岩塑性变形-混合型断裂-碎裂摩擦全过程；
S4、实验确定煤岩力学参数随超临界CO2作用时间的演化规律，输入材料参数并设置初始条件；
S5、设置增量步大小、输出参量、裂缝单元失效及接触对激活时刻的状态变量，分别通过有限元、离散元计算单元失效前后模型响应；
S6、分析超临界CO2作用不同时间下，受载煤岩试件变形-碎裂结果。
2.如权利要求1所述的一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法，其特征在于：所述步骤S1中，采用CT扫描方法确定煤岩中主要裂缝的分布，统计不同方位角裂隙出现的频数，统计相邻裂缝的平均间距及其方差；在此基础上编制主要裂缝自动生成程序，具体为：根据扫描确定的主要裂缝平均间距及其方差，在MATLAB中随机生成离散点，确定离散点的坐标，将散点连接成Delaunay三角形，由三角形相邻边的中垂线围成Voronoi多面体，生成多组Voronoi多面体，再优选出与CT扫描所得主要裂缝方位角结果接近的Voronoi多面体集合；由上述Voronoi多面体的各个面模拟煤岩中主要裂缝空间分布。
3.如权利要求1所述的一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法，其特征在于：所述步骤S2包括如下步骤：
S2.1、建立数值计算模型，划分实体单元网格；
在步骤S1的基础上，进一步根据数值模拟目的确定研究对象形状及尺寸，采用“切割”、“倒圆角”工具将Voronoi多面体模型截割成所需的数值模型；
网格划分前，依据下式确定单元最大尺寸L：
L＝9πEFe/[32(1-μ2)σp]
式中，E为弹性模量；μ为泊松比；Fe为断裂能；σp为抗拉强度；
通过Mesh模块“Seed edge”对数值模型布种，种子间距为L，进而通过“Seed part”命令对数值模型划分实体单元网格；
S2.2、在天然裂隙及实体单元交界处嵌入裂缝单元，并自动划分集合。
(1)更新实体单元节点编号；
网格划分后，通过CAE生成“.inp文件”，将文件内关键字“*nset”和“*elset”下方、“*asssembly”上方所有数字复制在Excel中，找到最大节点编号nmax、最大单元编号emax；
保持所有节点坐标不变，将第ni(0<i≤nmax)个节点复制a次，同时保持单元编号最小的节点编号ni,min不变，将其余单元的节点编号增加至最大节点编号的[(a-1)×10+ni,min]，其中，a为某个节点共用单元的个数；
(2)创建裂缝单元节点编号、单元编号；
若以四边形实体单元为例，记上述重新编号后的节点编号为N01,N02,N03,N04；共用第n0i个节点的相邻单元的编号为e01,e02,e03,e04。将新的节点编号和对应单元组合成有序数组：
Array＝(e01,e02,e03,e04,N01,N02,N03,N04)
将N01～N04逆序排列，即可得到四边形实体单元e01,e02,e03,e04之间裂缝单元的节点编号；
而裂缝单元最小的编号由实体单元最大编号所在数位+1得到的10的整数幂次方确定；
(3)在主要裂隙处嵌入裂缝单元；
在“.inp文件”中关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohelemAll”；按照第一列为裂缝单元编号、2～5列为节点号的顺序复制进相应的ei和Ni数字，从而在所有主要裂隙、所有实体单元边界处嵌入零厚度的裂缝单元；
(4)自动划分裂缝单元集合；
复制裂缝单元编号ei及其节点坐标Ni(xi,yi,zi)至Excel中，采用内积法求解所有裂缝单元的方向向量dcoh＝(d1,d2,d3)，以及Voronoi多面体各个面的方向向量dip，若二者共线，再通过空间向量方法计算裂缝单元上任意点到Voronoi面的距离，若为0，则判定裂缝单元位于Voronoi多面体的面上；将所有符合上述三个条件的裂缝单元编号、节点编号复制在.inp文件中，在关键字“*Part”之后“*Assembly”之前添加关键字
“*Elset＝cohElem01”，从而将所有主要裂隙上的裂缝单元组成一个名称为“cohelem01”的集合；
在此基础上，通过布尔运算获得所有实体单元边界处的裂缝单元集合“cohelem02”。
4.如权利要求1所述的一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法，其特征在于：所述步骤S3包括如下步骤：
S3.1、完整煤岩体弹塑性变形；
弹塑性变形阶段，各要素服从广义胡克定律：
σ＝E0:(ε-εp)
其中，σ为应力，E0为弹性刚度矩阵，ε为总应变。塑性应变
Figure FDA0003533011440000032
由屈服函数、塑性势函数确定，屈服函数F表达式为：
Figure FDA0003533011440000031
其中，F为屈服函数，<g>是Macaulay符号，
Figure FDA0003533011440000041
Figure FDA0003533011440000047
p＝trace(σ)/3，γ＝3(1-Kc)/(2Kc-1)，q＝[2(S:S)/3]1/2,S＝σ+pI式中，σi,max(i＝1,2,3)为最大有效主应力，
Figure FDA0003533011440000042
为双轴、单轴压缩屈服应力之比，σt与σc(εp)分别为拉压应力；Kc为控制屈服面在偏平面形状的参数，岩石类材料取值为0.667；
G＝[(δσttanψ)2+q2]-ptanψ
其中，δ为q-p平面内塑性势函数子午线在张拉区段的尖点曲率，一般取值为0.1；ψ为剪胀角；
S3.2、非贯通裂缝混合型断裂；
岩石断裂过程分为弹性变形和混合型断裂两个阶段；
弹性阶段本构关系为：
σc＝D0,cεc
一旦达到如下条件，进入韧性断裂过程：
Figure FDA0003533011440000043
其中，D0,c是弹性刚度矩阵，σc,n,σc,s,σc,t
Figure FDA0003533011440000044
是法向和两个切向方向的应力(峰值应力)，εc是应变；
为推导混合型韧性断裂本构方程，首先建立拉、剪、以及混合模式下断裂能Gc,n,Gc,S(Gc,S＝Gc,s+Gc,t，Gc,s与Gc,t为二个切向方向的断裂能),Gc,m表达式为：
Figure FDA0003533011440000045
其中，j表示n,S,m(即张拉、剪切、混合型)断裂模式；
断裂能混合比表达式为：
ξ＝Gc,S/(Gc,n+Gc,S)
煤岩完全断裂的判据为：
Figure FDA0003533011440000046
其中，
Figure FDA0003533011440000051
χ分别为完全断裂时的张拉、剪切型断裂能及材料参数；
拉、剪、混合型断裂模式下，各要素本构方程为：
Figure FDA0003533011440000058
其中，dc,j，Dc,j，σc,j，Sc,j分别为损伤变量、弹性模量、应力、塑性位移和总位移，且有：
Figure FDA0003533011440000052
其中，Sc,n与Sc,S分别为纯拉、剪条件下的位移；
假设混合型断裂模式下的塑性位移与剪切塑性位移相等，即有
Figure FDA0003533011440000053
其中，未知参数
Figure FDA0003533011440000054
通过以下步骤确定：①采用三点弯曲试验、加卸载直剪试验，获得拉、剪及混合模式下的σc,m-Sc,m曲线、σc,n-Sc,n曲线、σc,S-Sc,S曲线，计算相应的损伤演化曲线及断裂能
Figure FDA0003533011440000055
②通过“煤岩完全断裂判据表达式”对实验所得多组
Figure FDA0003533011440000056
数据进行拟合，得到相同材料参数χ下的断裂能混合比ξ；③在此基础上联立“峰值应力方程”，得到某个ξ下的峰值荷载
Figure FDA0003533011440000057
对应的σc,n与σc,S值，根据弹性表达式计算峰值位移分量Sc,n,Sc,S；④联立式“断裂能表达式”、“本构方程”，得到任意σc,m对应的σc,n与σc,S分量以及相对应的位移分量Sc,n,Sc,S；⑤将混合模式下的应力、位移替换为纯拉、纯剪的力学参数，并由“本构方程”计算σc,n，σc,S和Sc,n,Sc,S条件下的塑性应变、拉/剪损伤变量；⑦由“本构方程”计算得到得到任意混合比下的dc,m-Sc,m关系；由此确定混合型韧性断裂本构方程；
S3.3、贯通裂缝或块体的分离/挤压/剪切摩擦；
一旦断裂能满足断裂能表达式时，岩石完全断裂，在数值计算中删除该断裂面上的零厚度裂缝单元，消除其力学影响；
在此基础上，首先建立块体间分离、压缩、剪切摩擦关系判据：
①分离：相邻实体单元的任意二个节点间距l＞0，发生分离；
②挤压：l＝0且二者间存在压应力；
③剪切摩擦：若l＝0且二者间存在沿结构面的剪应力；
若相邻块体分离(l＞0)，二者间不存在相互作用力，块体运动服从牛顿第二定律；
若发生挤压，则结构面法向方向的力学本构关系为：
σn＝DnNmaxn/(Nmax-n)
其中，σn为压缩应力；Nmax为结构面最大闭合量，根据三维形貌扫描确定；n为结构面闭合量，变量；Dn为结构面法向模量，取值为其相邻岩块参数；
若发生剪切摩擦，粗糙结构面压剪摩擦本构方程：
Figure FDA0003533011440000061
其中，σS(σS,p)为(峰值)剪应力；DS为剪切刚度；S(sp)为(峰值)剪切位移；参数p由σS,p、Sp、残余剪应力σS,r和残余剪位移Sr确定；
步骤S3.4中，有限元FEM和离散元DEM耦合求解；
以被删除的零厚度裂缝单元为边界，将其内部分视作一个整体，根据部步骤S3.3方法，通过DEM方法求解所有离散块体的节点力；将此节点力作为边界条件，根据部步骤S3.1方法，通过FEM方法求解其内部完整块体的节点力；
通过实体单元与裂缝单元的共享节点传递力场数据，在此基础上，根据步骤S3.1计算实体单元的塑性应变、应力σ响应；根据步骤S3.2计算裂缝单元的损伤、位移、应力σc响应；若有断裂能
Figure FDA0003533011440000062
同时|σ-σc|≤10-5，则直接将实体单元、裂缝单元应力、位移响应传递至下一计算步骤；
若断裂能
Figure FDA0003533011440000063
表明裂缝单元完全断裂，此时删除零厚度裂缝单元，同时根据步骤S3.3计算煤岩块体间的分离/挤压/剪切摩擦；
若同时|σ-σc|＞10-5，则将附加节点力取值为|σ-σc|，根据步骤S3.1和S3.2重新计算实体单元、裂缝单元的节点力，直至|σ-σc|≤10-5成立；而后将实体单元、裂缝单元应力、位移响应传递至下一计算步骤，重复上述计算过程。
5.如权利要求1所述的一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法，其特征在于：所述步骤S4中，将煤岩试件在超临界二氧化碳(ScCO2)制取装置中浸泡0～60d，而后针对完整煤岩试件、含V型缺口的圆柱形试件、含粗糙贯通裂缝试件，分别采用“三轴压力试验机”、“电子万能试验机”、“多功能岩石三轴动态剪切试验机”，获取不同ScCO2浸泡时间下各试件的力学参数；
基于以上实验结果，在“Property”模块下设置1个场变量，场变量数值variables＝试件在ScCO2中的浸泡天数，同时设置相对应的“弹性模量”、“泊松比”、“拉/剪峰值应力”、“摩擦位移”；根据实验条件或工程条件，在“Load”模块下设置相应的应力边界、位移边界、初始应力场。
6.如权利要求1所述的一种超临界CO2作用下准脆性材料变形-碎裂的模拟方法，其特征在于：所述步骤S5中，在“Mesh”模块的“verify mesh”中，检查实体单元的最小稳定增量步，由此确定离散元DEM的最小增量时间步；
在“Mesh”模块下的Element type下设置裂缝单元失效及接触对激活时的变量；
在“Step”模块下的create field output中设置需要输出的变量，包括应力、应变、塑性应变、损伤变量、断裂能、摩擦应力、摩擦位移；
在此基础上，提交数值计算。
CN202210210428.3A 2022-03-04 2022-03-04 一种超临界co2作用下准脆性材料变形-碎裂的模拟方法 Pending CN114444230A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210210428.3A CN114444230A (zh) 2022-03-04 2022-03-04 一种超临界co2作用下准脆性材料变形-碎裂的模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210210428.3A CN114444230A (zh) 2022-03-04 2022-03-04 一种超临界co2作用下准脆性材料变形-碎裂的模拟方法

Publications (1)

Publication Number Publication Date
CN114444230A true CN114444230A (zh) 2022-05-06

Family

ID=81360369

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210210428.3A Pending CN114444230A (zh) 2022-03-04 2022-03-04 一种超临界co2作用下准脆性材料变形-碎裂的模拟方法

Country Status (1)

Country Link
CN (1) CN114444230A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117610395A (zh) * 2024-01-24 2024-02-27 西安交通大学 结晶岩压缩硬化记忆效应表征方法、装置、设备及介质

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117610395A (zh) * 2024-01-24 2024-02-27 西安交通大学 结晶岩压缩硬化记忆效应表征方法、装置、设备及介质
CN117610395B (zh) * 2024-01-24 2024-04-16 西安交通大学 结晶岩压缩硬化记忆效应表征方法、装置、设备及介质

Similar Documents

Publication Publication Date Title
Wang et al. Analyses and predictions of rock cuttabilities under different confining stresses and rock properties based on rock indentation tests by conical pick
Potyondy et al. Modeling of shock-and gas-driven fractures induced by a blast using bonded assemblies of spherical particles
Preece et al. A study of detonation timing and fragmentation using 3-D finite element techniques and a damage constitutive model
Hu et al. A distinct element based two-stage-structural model for investigation of the development process and failure mechanism of strainburst
CN113569442B (zh) 一种基于rkpm-pd耦合算法的岩石裂纹扩展预测方法
CN114444230A (zh) 一种超临界co2作用下准脆性材料变形-碎裂的模拟方法
CN114462124A (zh) 一种混凝土三维多相细观模型的建立与数值模拟方法
CN115859714B (zh) 一种基于fem-dem联合仿真的岩石爆破全过程模拟方法
da Silva et al. Analysis of the failure process by using the Lattice Discrete Element Method in the Abaqus environment
CN113836776A (zh) 一种爆破损伤预测模型构建方法
Gharehdash et al. Numerical investigation on fracturing of rock under blast using coupled finite element method and smoothed particle hydrodynamics
Bratov Incubation time fracture criterion for FEM simulations
Yahaghi et al. Development of a three-dimensional grain-based combined finite-discrete element method to model the failure process of fine-grained sandstones
Malhotra et al. An investigation on the accuracy of numerical simulations for single point incremental forming with continuum elements
CN115828472B (zh) 一种用于模拟滚磨光整加工工件表面残余应力的方法
Mirzaei Finite element analysis of deformation and fracture of cylindrical tubes under internal moving pressures
Zhao et al. Influence of inclination angles and confining pressures on mechanical behavior of rock materials containing a preexisting crack
Aliabadian et al. Simulation of dynamic fracturing of continuum rock in open pit mining
Ma Numerical modelling of underwater structural impact damage problems based on the material point method
Baranowski et al. Study of rock fracture under blast loading
Kurguzov Comparative analysis of failure criteria in building materials and rocks
Hao et al. A feasible approach for engineering-scale 3D blasting numerical modelling incorporating explosive charges and layout design
CN114818233A (zh) 一种煤层底板突水识别方法
Zidan et al. Modelling of damage patterns of RC concrete columns under demolition by blasting
Manko et al. Fracture mechanics of fractured rock masses and verification of rheological calculation models

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination