应用于岩石材料I型裂纹起裂和扩展过程微观分析方法与流程

本发明涉及一种裂纹起裂和扩展过程微观分析方法，尤其涉及一种应用于岩石材料i型裂纹起裂和扩展过程微观分析方法。

背景技术：

岩石的破坏过程是岩石工程安全控制的核心技术，岩石的微裂纹扩展机理和控制成为岩石工程研究的前沿热点问题。由于岩石的微裂纹起裂机理不清楚，给岩石工程工作者带来很大困惑。岩石是多孔介质，但其骨架成分是固体材料，很多都是结晶体构成，其中二氧化硅是一种典型的分子组成成分，研究晶体的断裂过程是研究微裂纹起裂和扩展的关键技术之一。也是更好的理解高应力、高温度等复杂环境下岩石材料脆裂蠕变的基础。目前中国还没有采用分子动力学理论进行岩石材料起裂和扩展机理分析控制的研究成果，因此基于分子动力学理论建立岩石材料裂纹起裂和扩展的微观分析方法具有重要意义，也是寻求控制微裂纹扩展方法的必经之路，具有重要科学意义和实用价值。

技术实现要素：

发明目的：本发明的目的在于提供一种应用于岩石材料受拉条件下微裂纹起裂和扩展控制的分析方法。

技术方案：一种应用于岩石材料i型裂纹起裂和扩展过程微观分析方法，该方法包括以下步骤：

(1)采用lammp软件建立sio2裂纹起裂的微观模型；

(2)设置初始模型温度、初始裂纹长度及位移加载速度，模拟分析裂纹扩展过程，并在裂纹扩展路径上取贴近断裂面的相同间距的点，点的个数≥6；根据模拟结果，计算裂纹一端区域的平均应力及裂纹总长度；

(3)根据步骤(2)中数据计算裂纹应力强度因子，裂纹应力强度因子按照断裂力学应力强度因子公式来计算：

其中，ki是裂纹应力强度因子，y是与裂纹形状和加载方式有关的系数，这里取1.12，σ是a区域平均应力，c是裂纹总长度；

(4)根据步骤(3)中数据计算起裂裂纹驱动力，裂纹驱动力的表达式为：

其中，gi是裂纹驱动力，b是模型宽度，是指系综总势能随着裂纹扩展降低的量，e是弹性模量；

(5)建立微观裂纹扩展平均速率的计算方法：在裂纹扩展的方向上选择直径为6倍以上分子长度的圆，圆的个数≥3，通过软件读取每个圆的能量历时变化过程，以及各个点的位移历时变化规律，根据选择圆的能量变化突变计算裂纹扩展平均速率，计算裂纹扩展速率与能量之间的关系；

(6)分析不同系综、不同温度、不同加载速率、不同力加载或不同裂纹板厚度对裂纹扩展过程的影响。

其中，步骤(1)优选采用vmd建模软件建立模型后导入lammp软件平台。所述sio2的模拟采用vashishta势函数，表达式为：

式中，表示两体相互作用，表示三体相互作用，r表示原子之间的距离，m，n是指定的原子下标，z表示原子电荷，α表示原子的极化性，σ表示最接近点的距离，θ是海维赛德函数。

步骤(2)中所述初始模型温度为10k，初始裂纹长度位移加载速度为所述相同间距的点的个数优选为11个；所述应力的计算采用维里应力，维里应力表征原子间相互作用力，可正可负，表达式为：

其中，v是需要计算应力部分的体积，i，j表示原子i周围的j原子，表示i原子和j原子在α方向的位置，表示β方向j原子对i原子的力，mⁱ表示i原子的质量，和分别表示i原子在α和β方向的速度。

步骤(5)中所述圆的个数优选为3个。

有益效果：与现有技术相比，本发明的优点是：(1)首次提出采用分子动力学模拟岩石材料的不同系综、不同加载方式、不同温度下微裂纹起裂和扩展的过程，建立了微观分子动力学模拟下应力强度因子的计算方法；概念清晰，模拟结果直观，容易被岩石工程学者接受、理解和分析，为微观研究岩石材料的起裂和扩展机理提供有效新途径。(2)基于分子动力学模拟技术，从能量角度分析裂纹扩展的过程，更符合实际情况，弥补了断裂力学的宏观计算分析概念太抽象的缺陷，并建立了微观模拟中裂纹扩展速率的计算方法，为深入分析岩石脆裂蠕变等微裂纹扩展全过程及其机理研究提供了理论基础。

附图说明

图1为sio2分子的空间构型；

图2为sio2计算模型示意图；

图3为裂纹扩展过程示意图；

图4为裂纹扩展速率与长度曲线；

图5为a区域应力与裂纹应力强度因子曲线；

图6为微裂纹扩展过程的能量变化曲线；

图7为弹性力学裂纹尖端应力场示意图；

图8为裂纹扩展过程裂尖应力的变化曲线；

图9为nvt系综下裂纹扩展过程的计算分析结果；

图10为温度对裂纹扩展的影响规律对比分析曲线；

图11为加载速率对裂纹扩展的分析曲线；

图12为力加载模式下裂纹扩展过程分析曲线；

图13为裂纹板厚度对裂纹扩展过程分析曲线。

具体实施方式

下面结合附图对本发明的技术方案作进一步说明。

一种应用于岩石材料i型裂纹起裂和扩展过程微观分析方法，基于二氧化硅建立了微裂纹起裂和扩展过程的分子动力学数值模拟，分析了不同系综、不同温度、不同加载速率、不同力加载或不同裂纹板厚度对裂纹扩展过程的影响规律。在lammp开放平台基础上编制了响应程序。形成了一种基于分子动力学模拟技术的岩石材料受拉条件下微裂纹扩展过程分析的新方法。

本发明的具体分析过程如下：

针对岩石基材二氧化硅(sio2)基于lammp开源软件平台，编制程序建立微裂纹扩展分析分子动力学模拟模型，进行不同工况下裂纹扩展规律分析，为微裂纹的控制提供理论基础。

1、势函数与模型

sio2的模拟采用vashishta势函数，表达式为：

式中，表示两体相互作用，表示三体相互作用，r表示原子之间的距离，m，n是指定的原子下标，z表示原子电荷，α表示原子的极化性，σ表示最接近点的距离，θ是海维赛德函数。

根据实验获得vashishta势函数的参数如表1-4所示：

表1si-si-si参数表

表2o-o-o参数表

表3o-si-si参数表

表4si-o-o参数表

二氧化硅(sio2)的分子动力学模拟模型采用vmd软件建立，具体分子空间构型如图1所示，其中图1(a)为单个sio2分子的空间构型，图1(b)为多个sio2分子的空间构型，模型实例如图2所示，尺寸为

2微裂纹扩展过程及其力学分析

2.1微裂纹扩展过程模拟结果分析

设置初始模型温度为10k，初始裂纹长度位移加载速度为模拟分析裂纹扩展过程如图3所示。模拟结果显示：拉伸过程中，能量逐渐向裂纹尖端聚集，裂纹尖端出现应力集中，势能较小的原子之间慢慢被拉开，直至势函数作用距离之外，尖端开裂。如图3(a)所示，在加载10ps左右，裂纹尖端出现空位，显然，该空位没有出现在裂纹正前方，是二氧化硅的构型决定的，即原子的空间排布导致空位处的原子之间势能较低。加载到20ps，如图3(b)所示，预制裂纹尖端及空位各自形成一个裂纹尖端，即“裂纹1”和“裂纹2”，随后，在两个裂纹尖端的“竞争”中，裂尖之间的原子被拉开直至拉断，裂尖1被拉开，裂纹沿着裂尖1向前扩展，如图3(c)(d)所示，其中，图3(c)加载到31ps，图3(d)加载到70ps，之后加载到170ps，如图3(e)所示，图3(f)为加载到263ps。图3(d)(e)(f)分别表示裂纹扩展至图2(a)所示的“圆1”、“圆2”和“圆3”的位置。

在裂纹扩展路径上取11个贴近断裂面的点，如图4(a)所示，这些点y方向位移变化如图4(b)所示，判断各点起裂时间，并计算裂纹长度以及裂纹扩展的大致速度变化，如图4(c)所示，由图4(c)可知裂纹在扩展过程中并不是匀速前进的，而是上下震荡形式的，结合图4(d)可看出，裂纹长度曲线斜率在260ps之前逐渐变小，之后曲率增大。所以速率变化的总体趋势是先减小(260ps之前)后增大。不计裂纹停止的时间，裂纹扩展速率在(即200m/s)左右，最高速率超过(即400m/s)，该计算结果与s.chandra等人计算的al在1k温度下扩展速率320m/s比较接近，所以本模型计算结果基本合理可靠。

加载过程中，如图4(d)所示，整个过程应力变化为先减小后略微增大，从260ps左右开始再次减小。应力的计算采用的是维里应力，维里应力表征原子间相互作用力，可正可负，表达式为：

其中，v是需要计算应力部分的体积，i，j表示原子i周围的j原子，表示i原子和j原子在α方向的位置，表示β方向j原子对i原子的力，mⁱ表示i原子的质量，和分别表示i原子在α和β方向的速度。

2.2微裂纹的应力强度因子计算方法

初始应力较高的原因是，裂纹扩展前，a区域的右下角处于预制裂纹尖端的附近，该处应力值比较大，导致a区域总体平均应力较大。裂纹扩展后，裂纹尖端逐渐远离a区域，该区域应力下降，应力峰值时刻对应裂纹加速扩展时刻，裂纹应力强度因子按照断裂力学应力强度因子公式来计算：

其中，ki是裂纹应力强度因子，y是与裂纹形状和加载方式有关的系数，这里取1.12，σ是a区域平均应力，取图5(a)中的数据，c是裂纹总长度，取图4(d)中的数据。应力强度因子计算结果如图5(b)所示，应力强度因子逐渐升高，20ps左右裂纹起裂，此时裂纹应力强度因子可能达到了断裂韧性导致裂纹开裂。裂纹驱动力的表达式为：

其中，gi是裂纹驱动力，b是模型宽度，是指系综总势能随着裂纹扩展降低的量，e是弹性模量。20ps左右裂纹起裂，二氧化硅弹性模量约为70gpa，计算得到起裂裂纹驱动力为0.0057j/m²。

2.3建立微观裂纹扩展平均速率的计算方法

如图6(a)所示，在裂纹扩展的方向上选择直径为6倍以上分子长度的3个圆，图6(b)为指定圆内的能量，图6(c)为模型系统的能量。读取图6(b)中能量增加过程中的突变点a、b、c对应的时间，利用原中心点的距离，可计算出裂纹从圆1扩展到圆2的平均速率以及从圆2扩展到圆3的平均速率

裂纹扩展过程中，图6(a)中三个圆形裂纹尖端区域能量变化如图6(b)所示，三个区域变化趋势大致相同，在裂纹扩展的瞬间能量突然升高，随后维持稳定，裂开前后动能变化不大，而且nve系综为封闭系统，不与外界发生能量交换，总能量主要包括动能、势能和表面能，不考虑其他能量耗散项，因此可认为能量升高是因为三个区域都包含裂纹断裂后的表面，储存有表面能。从图6(c)中可看出，系统动能变化不明显，可以认为三个圆形裂纹尖端区域能量升高的部分均是表面能，因此可以计算出表面能，每个圆形区域裂纹开裂截面面积为开裂后能量增长约为0.6ev，所以表面能为0.23j/m²。计算过程是在极低温度(10k)下进行的，而且断裂过程裂纹尖端始终是尖锐的，所以模型是高度脆性的，根据lawn给出的高度脆性材料的表面能在10^-1～10范围内，本模型计算的表面能结果在该范围内。

2.4微裂纹起裂机理分析

根据断裂力学理论，裂纹尖端应力场如图7所示，i型张拉裂纹尖端应力和位移的弹性平面应力解析表达式为：

式中，σ1、σ2分别是裂纹尖端附近某一点的主应力，r是距离尖端原点的距离，θ是与x轴方向的夹角。r0处的裂纹张开量为δ，一般认为上述公式主要适用于裂纹尖端主要区域r1<r<r2内，但区域的上下边界仍存在争议。

由上述公式可知，距离裂纹尖端越近，应力越大，理论上，裂纹尖端应力趋于无穷大，这是不符合物理规律的，所以裂纹尖端应力到底多大时裂纹开裂是弹性力学解无法触及的“黑匣子”。物理层面，裂纹尖端是由原子分布形成的，图3说明了裂纹扩展是原子之间相互作用断开即键断裂导致的，所以裂纹尖端应力达到一定程度，原子键被拉断，裂纹扩展。为了探究裂纹尖端应力变化，得到裂纹扩展过程，计算图6(a)中圆1、2、3三个区域水平和竖直方向应力变化，如图8所示。微观尺度上，原子时刻都在运动，因此应力值是波动的。图8(a)是三个裂纹尖端区域水平方向应力变化，裂纹依次进入三个区域，三个区域先后出现峰值，峰值对应裂纹进入该区域的时间，可见裂纹扩展过程尖端水平向应力最大可达到4gpa，同理，8(b)表示裂纹扩展过程中三个尖端区域竖直向应力变化，裂纹进入的瞬间应力可高达6gpa左右。对比(a)和(b)可发现，水平向应力增加是逐渐增加，逐渐降低，垂直向则是快速增加，裂纹扩展后几乎是瞬间下降至裂纹未靠近的应力值范围，说明裂纹尖端向前推移的过程中，水平向应力受影响的原子远多于垂直向应力受影响的原子。裂纹扩展时，裂纹尖端应力的“奇异性”是不存在的，应力增加到6gpa左右裂纹会裂开，向前扩展。

3不同条件对裂纹扩展行为的影响

3.1系综对裂纹扩展行为的影响

在nvt系综下进行计算，计算结果如图9所示，图9(a)为nve和nvt系综总能量变化，由图9(a)可知nvt总能量增长幅度小于nve系综，这与nvt系综可以与外界环境有能量交换导致的。图9(b)为nve和nvt系综a区域应力变化，图9(b)可知两个系综a区域y方向应力开始时刻相似，随后nvt系综下应力增大，nve系综基本不变，最后二者又都经历减小的过程。图9(c)为nve和nvt系综裂纹长度变化，图9(d)为nve和nvt系综裂纹扩展速率，由图9(c)、(d)可知nve系综下裂纹扩展平均速率高于nvt系综。图9(e)为nve和nvt系综裂纹应力强度因子，由图9(e)可知nvt系综下，裂纹扩展的应力强度因子大于nve系综，这是由于nvt系综允许与外界能量交换，导致裂纹扩展需要更大的裂纹应力强度因子。

3.2不同温度对裂纹扩展行为的影响

在nve系综，基本条件不变，按位移速度加载，分别计算10k、20k和30k温度下裂纹扩展行为如图10所示，其中图10(a)为不同温度系统总能量；图10(b)为不同温度a区域应力；图10(c)为不同温度圆1区域能量变化；图10(d)为不同温度圆2区域能量变化；图10(e)为不同温度圆3区域能量变化。由图10(a)(b)可看出，系综总能量在扩展后期，温度较高的模型由于动能也稍高，所以总能量偏高，a区域平均应力变化基本无变化。对比图10(c)(d)和(e)，可发现10k，20k和30k温度下裂纹扩展路径上三个圆区域能量变化基本一致。

3.3不同加载速率对裂纹扩展行为的影响

在nve系综、10k温度下，将加载速率由降低为其裂纹扩展行为如图11所示，其中，图11(a)为空位形成示意图；图11(b)为裂纹沿裂纹尖端扩展示意图；图11(c)为不同加载速率圆1能量变化；图11(d)为不同加载速率圆2能量变化；图11(e)为不同加载速率圆3能量变化；图11(f)为不同加载速率系统总能量变化；图11(g)为a区域平均应力变化；图11(h)为系统温度变化。由图11(a)(b)可看出，初始裂纹尖端行为与图3(a)相似，但是随后低加载速率下裂纹并不沿着空位扩展，而是从裂纹尖端向前扩展。原因是虽然空位处原子间结合能较低，但是慢速加载导致原子有足够时间重新调整位置，降低势能，提高整体稳定性，空位得到填充，于是裂纹沿着尖端扩展。由图11(c)(d)(e)可知，圆2区域在裂纹扩展后能量变化稍低于圆1和3区域，裂纹穿过圆2区域用时也比较长，该区域原子有更多的时间调整位置，降低势能，提高稳定性，所以最后总能量偏低。总体势能、动能变化，a区域平均应力变化，以及总体温度变化与加载速率计算结果规律基本相同，但加载速率可能会引发裂纹扩展路径的变化，这与实际试验结果基本一致。

3.4不同力加载对裂纹扩展行为的影响

在nve系综、10k温度下，将图3(a)中速度加载换成力加载，并将圆2、3位置调整一下，调整之后如图12所示，其中图12(a)为力加载模型示意图；其中图12(b)为裂纹沿裂纹尖端扩展(5ps)示意图；其中图12(c)为纵向起裂(25ps)示意图；其中图12(d)为纵向裂纹扩展(26ps)示意图；其中图12(e)为圆1水平应力与竖直应力；其中图12(f)为圆2水平应力与竖直应力；其中图12(g)为圆3水平应力与竖直应力；其中图12(h)为温度变化曲线；其中图12(i)为动能、势能及总能量变化曲线；其中图12(j)为三个圆区域能量变化曲线。调整加载力的大小，当加载应力小于1.793gpa，裂纹不扩展，加载应力大于等于1.793gpa，裂纹扩展，初始裂纹长度为根据如下公式可计算出起裂应力强度因子与裂纹扩展驱动力分别为：

比速度加载情况下计算的裂纹扩展应力强度因子大4.9倍，裂纹驱动力增大24倍。裂纹一旦开裂，迅速扩展，由于裂纹扩展速度过快，而且是端部加载，导致模型纵向断裂，如图12(d)所示，究其原因，是圆2、3部分的水平应力大于竖直应力，如图12(f)(g)；断裂后，断裂的部分受力快速运动，所以系统动能急剧升高，系统温度也急剧升高，如图12(h)(i)所示。力加载方式显示，力小于某阈值时，微裂纹不会扩展，这与断裂力学理论一致。

3.5不同模型厚度对裂纹扩展行为的影响

在nve系综、10k温度下，初始裂纹长度为加载速度为模型厚度由原来的增大到裂纹板厚度对裂纹扩展过程的影响如图13所示，其中图13(a)为模型正视图；图13(b)为模型左视图；图13(c)为系统能量变化曲线；图13(d)为a区域应力变化曲线；图13(e)为不同厚度圆1能量变化曲线；图13(f)为不同厚度圆2能量变化曲线；图13(g)为不同厚度圆3能量变化曲线；图13(h)为不同厚度应力强度因子。研究结果发现，模型厚度会改变裂纹扩展行为。如图13(a)(b)，可以看出在扩展的最后阶段，裂纹尖端偏离了原来的方向，有一层原子“粘”在下表面，裂纹沿着新的表面扩展，分析认为产生这一现象的原因主要是裂纹尖端始终自洽寻找断裂韧度或者扩展阻力小的地方扩展，厚度方向的增加使得裂纹扩展过程中裂纹尖端的“弱键”增加，使得扩展行为发生了细微的变化。图13(d)可看出，厚度增加后，裂纹扩展过程a区域平均拉应力增加，曲线形状更像nvt系综下的应力变化曲线。厚度增加后，系统内原子数量增加，所以系统动能势能和总能量变化均增大，如图13(c)所示，裂纹扩展形成的表面增大，扩展后三个圆区域内能量增幅变大，如图13(e)(f)(g)所示。模型厚度增大导致开裂过程中所需应力强度因子增大，如图13(h)所示，这与断裂力学实验结果基本一致，但分子动力学模拟分析方法更清楚的揭示了断裂裂纹边界条件对裂纹扩展的影响机理。

由此可见，采用本发明的分析方法可以更全面的揭示岩石材料微裂纹起裂和扩展的影响机理，对于不同基质的岩石材料采用同样的方法，同样可以更全面的揭示微裂纹扩展机理，从而为控制微裂纹的扩展提供了理论依据。