一种基于分子动力学模拟求解矿物微观弹性参数的方法

：
1.本发明涉及分子动力学模拟技术，具体涉及的是一种基于分子动力学模拟求解矿物微观弹性参数的方法。


背景技术：

2.随着常规油气的大规模勘探开发，其增储上产难度也越来越大，这与人们不断增长的油气资源需求之间的矛盾日益突出。我国页岩油气资源储量相当丰富，分布广泛，发展潜力巨大，若能合理利用页岩油气资源可很大程度上解决我国油气能源问题。页岩储层与常规的砂岩和碳酸盐储层有所不同，页岩储层既是烃源岩也是储集岩，且基质孔隙度小，孔喉半径小，渗透率极低，自然条件下基本没有工业产能，不具备商业开采价值，形成页岩油气工业产能的主要手段是对页岩储层进行大规模压裂改造，压裂效果受到储层岩石物理力学特性及储层地质条件等多方面影响，因此开展页岩力学性质研究意义重大。岩石力学参数主要包括体积模量、剪切模量、杨氏模量和泊松比等，这些力学参数可通过实验得出，例如根据不同围压下三轴抗压实验能够确定岩样的内摩擦角和内聚力；利用岩石单轴抗压实验可以绘制应力应变曲线，从而能够计算弹性模量、泊松比及抗压强度等力学参数。
3.页岩由多种矿物组成，具有极强的非均质性、各向异性和多尺度属性，现有的室内试验及理论分析方法确定岩石力学参数仍有较大局限性，目前岩石力学参数的求取方法主要为实验室岩心测定法和测井资料解释方法，精度较高的岩心测定方法不足之处在于测定过程复杂、消耗成本较大；从岩石的声学特性计算岩石力学参数虽然成本较低，且同时可得到井点处连续的岩石力学参数曲线，但只能进行一维点上的参数预测而无法进行空间预测。宏观力学性质本质上由其微观结构和各矿物成分力学性质决定，因此可以从微结构的本质属性出发解释宏观的现象，针对上述问题，分子动力学模拟方法是一种很好的选择。
4.分子动力学基于经典力学，利用牛顿运动定律模拟体系内原子或分子的运动轨迹，并通过对其不同状态所构成的系综进行统计平均计算体系的结构和性质。这种方法能够呈现不同物质之间相互作用的微观细节且能够准确复现实验结果，因此在物理化学、生命科学和材料等领域均有长远的发展，石油相关领域也逐渐开始应用这种技术，例如通过平衡分子动力学方法分析纳米孔隙内烷烃与水的结构特征、扩散系数等，以此明确油水两相赋存状态；通过非平衡分子动力学方法分析纳米孔隙内烷烃与水的速度、滑移、及流量，以此明确油水两相流动规律等相关研究。目前利用分子动力学方法研究纳米尺度下物质的力学性能和变形机理已有一定进展，例如通过仿真碳纳米管计算其力学性能参数，如应力应变关系、弹性模量等；建立弯曲加载下α-al2o3纳米线模型，计算不同加载速率下α-al2o3纳米线的原子应力和应变，揭示了加载速率对其弯曲力学行为的影响规律等研究。本发明利用分子动力学模拟方法构建石英、伊利石、方解石矿物模型，选取clayff以及xiao等人开发的力场参数并予以验证后求解了各页岩主要矿物的弹性常数及相关力学参数。


技术实现要素：

5.本发明的目的是一种基于分子动力学模拟求解矿物微观弹性参数的方法，这种基于分子动力学模拟求解矿物微观弹性参数的方法用于解决现有的室内试验及理论分析方法确定岩石力学参数仍有较大局限性的问题。
6.本发明解决其技术问题所采用的技术方案是：这种基于分子动力学模拟求解矿物微观弹性参数的方法包括以下步骤：
7.从微观尺度分别建立石英、伊利石及方解石的分子动力学模型，使用msi2lmp程序将各分子动力学模型坐标转换为lammps识别的数据文件；
8.分别描述石英、伊利石及方解石内部相互作用，通过对比驰豫前后矿物成分的晶格参数，确定所选力场合理性；
9.通过对分子动力学模型不同方向施加形变，收集每次应变后应力应变数据，准确求解石英、伊利石及方解石的弹性常数、体积模量、剪切模量；采用delta方式对分子动力学模型模型施加应变，定义形变量为1.0
×
10-6
与施加载荷方向晶格长度之积，原子随机抖动量为分别对x、y、z、xy、xz、yz六个方向施加拉伸载荷和压缩载荷，每次操作后使用能量最小化的方式松弛原子位置以确保模型结构稳定，其中能量的停止容差为1.0
×
10-8
，力的停止容差为1.0
×
10-10
，最大迭代次数为100，计算力或能量的最大次数为1000；不同方向进行拉伸或压缩操作后获得新的应力张量，共进行12次计算求解全部弹性常数，进一步计算石英、伊利石及方解石的体积模量、剪切模量。
10.上述基于分子动力学模拟求解矿物微观弹性参数的方法具体为：
11.1)分别建立石英、伊利石、方解石分子模型，利用msilmp程序将各模型坐标转换为lammps识别的数据文件；
12.2)描述石英和伊利石原子间相互作用，描述方解石原子间相互作用，设置非键势能、键伸缩能、键角弯曲能以及不同原子电荷参数，各公式如下：
13.总相互作用能e：
14.e＝e
nobonded
+e
bonds
+e
angles
15.非键结势能e
nobonded
(包括l-j势能和静电势能)：
[0016][0017]
键伸缩势能e
bonds
：
[0018][0019]
键角弯曲势能e
angles
：
[0020][0021]
式中：ε
ij
、σ
ij
为原子i和j之间的能量参数和距离参数；r
ij
为分离距离；kr和k
θ
为力常数；θ为键的夹角；e为电子电荷；qi为i原子的部分电荷；qj为j原子的部分电荷；∈0为真空介电常数(8.85419
×
10-12
f/m)；r和θ分别为键长和键角；r0和θ0为相应量的平衡值；εi和σi为常见的lennard-jones能量和尺寸参数；
[0022]
3)设置初始模拟参数，计量单位选择real格式，截断半径为边界条件为三维周期性边界，计算精度设为1
×
10-4
，时间步长为1fs；
[0023]
4)使用共轭梯度分别对石英、伊利石、方解石分子模型进行能量最小化，对模型初始构型进行调整；
[0024]
5)定义热力学输出信息，包括步数、温度、压力、总势能、动能，晶格长度、晶格角度，设置每间隔100步输出一次；
[0025]
6)采用nose hoover方法控温，使各模型先后在nvt、npt系综下充分驰豫，驰豫时间为50ps，设置温度为348.73k，压力为45.82mpa，在npt系综下输出轨迹文件，每100步输出一次；
[0026]
7)对比驰豫前后石英、伊利石及方解石各自的晶格常数及矿物结构稳定性，验证所选力场合理；
[0027]
8)对石英、伊利石及方解石结构单元执行扩胞操作，将扩胞至合适尺寸的矿物模型转换为数据文件；
[0028]
9)定义形变量为施加应变方向晶格长度
×
10-6
，原子随机抖动量为选取验证合理后的相关势函数；
[0029]
10)使用能量最小化的方式松弛原子位置以确保模型结构稳定，其中能量的停止容差为1.0
×
10-8
，力的停止公差为1.0
×
10-10
，最大迭代次数为100，计算力或能量的最大次数为1000，计算模型初始长度及变化长度，以此来定义应变；
[0030]
11)使用delta方式分别对各模型x、y、z、xy、xz、yz六个方向施加拉伸载荷和压缩载荷，每次仅取施加载荷方向的应变，共进行12次数据收集；
[0031]
12)将拉伸过程得到的弹性常数与压缩过程得到的弹性常数结果做平均处理，得到弹性常数；
[0032]
13)根据上述求得的弹性常数进一步求解矿物的体积模量、剪切模量，石英和方解石为三角晶系，伊利石为单斜晶系，三角晶系和单斜晶系对应的弹性矩阵和力学参数计算公式如下:
[0033]
三角晶系具有6个独立弹性常数，弹性范围内，其应力σ与应变ε关系通过广义胡克定律表示：
[0034][0035]
相应的力学参数求解公式如下：
[0036]
体积模量b：
[0037]
b＝(2c
11
+c
33
+2c
12
+4c
13
)/9
[0038]
剪切模量g：
[0039]
g＝(7c
11
+2c
33
+12c
44-5c
12-4c
13
)/30
[0040]
杨氏模量e：
[0041][0042]
泊松比v：
[0043][0044]
单斜晶系具有6个独立弹性常数，其应力应变关系表示为：
[0045][0046]
式中：σ
xx
、σ
yy
、σ
zz
为x、y、z三方向主应力；τ
yz
、τ
zx
、τ
xy
为yz、zx、xy三方向切应力；ε
xx
、ε
yy
、ε
zz
为x、y、z三方向主应变；γ
yz
、γ
zx
、γ
xy
为yz、zx、xy三方向切应变；c
11-c
66
为弹性矩阵元。相应的力学参数求解公式如下：
[0047]
体积模量b：
[0048][0049]
剪切模量g：
[0050][0051]
单斜晶系中杨氏模量与泊松比的计算公式与三角晶系中杨氏模量与泊松比的计算公式相同。
[0052]
有益效果：
[0053]
1.本发获得的矿物微观弹性参数，结合等效夹杂理论与有限元分析方法即可获得具有各种矿物含量的不同岩性岩石的宏观力学性质，为页岩力学性质的求取提供了一种更为简便的方法，也为目前页岩取芯困难、标准岩样制备成功率低等导致页岩力学性质测定困难等问题提供了有效解决手段。
[0054]
2.本发明从弹性矩阵元入手求解矿物力学参数，解决了各向异性矿物不同方向力学性质不同的问题；计算弹性常数时综合考虑各个方向力学性质，相比施加单一方向压缩或拉伸所得到的数据更为精确。。
[0055]
3.本发明通过对比驰豫前后各矿物分子模型晶格参数及模型结构稳定性来验证力场，证明了应用的势参数合理可靠，有助于验证分子动力学模拟准确性，对提高分子动力学模拟结果的可信度以及分子动力学技术的广泛应用具有重要意义。
[0056]
4.本发明使用分子动力学模拟方法研究页岩力学性质，一方面可以节约实验成本；另一方面宏观力学性质本质上由其微观结构和各矿物成分力学性质决定，该方法基于分子尺度对页岩微结构施加载荷，分析应力应变并计算力学参数，为页岩宏观力学性质研
究奠定理论基础。
附图说明：
[0057]
图1为页岩主要矿物成分分子结构模型图，图中(a)石英；(b)伊利石；(c)方解石；
[0058]
图2为扩胞后模型结构图，图中(a)石英；(b)伊利石；(c)方解石；
[0059]
图3为分子动力学模拟矿物力学性质流程图。
具体实施方式：
[0060]
下面结合附图对本发明做进一步说明：
[0061]
这种基于分子动力学模拟求解页岩矿物微观力学参数的方法，具体为从微观尺度建立主要矿物成分(石英、伊利石、方解石)的分子动力学模型，选择clayff力场参数描述石英、伊利石内部相互作用，选取xiao等人开发的力场参数描述方解石内部相互作用，通过对比驰豫前后矿物成分的晶格参数，确定所选力场合理性。设定合理模拟参数及方法步骤，通过对模型不同方向施加形变，收集每次应变后应力应变数据，准确求解页岩主要矿物成分的弹性常数及相关力学参数。主要内容如下：分别构建石英、伊利石、方解石分子结构模型，大规模原子/分子并行模拟器能够识别数据文件，故使用msi2lmp程序将矿物模型坐标进行转换，编写程序从而实现对矿物相关性能的分子动力学模拟计算，力场参数如下：
[0062]
石英、伊利石势参数选取
[0063][0064]
方解石势参数选取
[0065][0066]
力学求解参数部分采用delta方式对模型施加应变，定义形变量为1.0
×
10-6
与施加载荷方向晶格长度之积，原子随机抖动量为分别对x、y、z、xy、xz、yz六个方向施加拉伸载荷和压缩载荷，每次操作后使用能量最小化的方式松弛原子位置以确保模型结构稳定，其中能量的停止容差为1.0
×
10-8
，力的停止容差为1.0
×
10-10
，最大迭代次数为100，计算力或能量的最大次数为1000。不同方向进行拉伸或压缩操作后获得新的应力张量，共进行12次计算求解全部弹性常数，从而进一步计算不同矿物体积模量、剪切模量等相关力学参数。
[0067]
更具体的，这种基于分子动力学模拟求解矿物微观弹性参数的方法包括如下步骤：
[0068]
1)建立主要矿物分子模型(石英、伊利石、方解石分子结构见附图1)，利用msilmp程序将模型坐标转换为lammps识别的数据文件。
[0069]
2)选取势参数，使用clayff力场描述石英和伊利石原子间相互作用，使用xiao等人开发的力场描述方解石原子间相互作用，设置非键势能、键伸缩能、键角弯曲能以及不同原子电荷等参数，相关势能公式如下：
[0070]
总相互作用能：
[0071]
e＝e
nobonds
+e
bonds
+e
angles
[0072]
非键结势能(包括l-j势能和静电势能)：
[0073][0074]
键伸缩势能：
[0075][0076]
键角弯曲势能：
[0077][0078]
式中：ε
ij
、σ
ij
为原子i和j之间的能量参数和距离参数；r
ij
为分离距离；kr和k
θ
为力常数；θ为键的夹角；e为电子电荷；qi为i原子的部分电荷；qj为j原子的部分电荷；∈0为真空介电常数(8.85419
×
10-12
f/m)；r和θ分别为键长和键角；r0和θ0为相应量的平衡值；εi和σi为常见的lennard-jones能量和尺寸参数；
[0079]
3)设置初始模拟参数，计量单位选择real格式，截断半径为边界条件为三维周期性边界，计算精度设为1
×
10-4
，时间步长为1fs。
[0080]
4)使用共轭梯度对模型进行能量最小化，从而对模型初始构型进行调整，实现结构优化。
[0081]
5)定义热力学输出信息，包括步数、温度、压力、总势能、动能，晶格长度、晶格角度等，设置每间隔100步输出一次。
[0082]
6)采用nose hoover方法控温，使模型先后在nvt、npt系综下充分驰豫，驰豫时间为50ps，设置温度为348.73k，压力为45.82mpa，在npt系综下输出轨迹文件，每100步输出一次。
[0083]
7)对比驰豫前后各矿物的晶格常数及矿物结构稳定性，以验证所选力场合理(附表2)。
[0084]
8)对矿物结构单元执行扩胞操作(附图2)，将扩胞至合适尺寸的矿物模型转换为数据文件。
[0085]
9)定义形变量为施加应变方向晶格长度
×
10-6
，原子随机抖动量为选取验证合理后的相关势函数。
[0086]
10)使用能量最小化的方式松弛原子位置以确保模型结构稳定，其中能量的停止容差为1.0
×
10-8
，力的停止公差为1.0
×
10-10
，最大迭代次数为100，计算力或能量的最大次数为1000，计算模型初始长度及变化长度，以此来定义应变。
[0087]
11)使用delta方式分别对模型x、y、z、xy、xz、yz六个方向施加拉伸载荷和压缩载荷，每次仅取施加载荷方向的应变，共进行12次数据收集。
[0088]
12)将拉伸过程得到的弹性常数与压缩过程得到的弹性常数结果做平均处理。
[0089]
13)根据上述求得的弹性常数可进一步求解矿物的体积模量、剪切模量等力学参数，求三种矿物力学参数验证结果和矿物力学参数汇总。石英和方解石为三角晶系，伊利石为单斜晶系，三角晶系和单斜晶系对应的弹性矩阵和力学参数计算公式如下:
[0090]
三角晶系具有6个独立弹性常数，弹性范围内，其应力σ与应变ε关系可通过广义胡克定律表示：
[0091][0092]
相应的力学参数求解公式如下：
[0093]
体积模量b：
[0094]
b＝(2c
11
+c
33
+2c
12
+4c
13
)/9
[0095]
剪切模量g：
[0096]
g＝(7c
11
+2c
33
+12c
44-5c
12-4c
13
)/30
[0097]
杨氏模量e：
[0098]
[0099]
泊松比v：
[0100][0101]
单斜晶系具有6个独立弹性常数，其应力应变关系可表示为：
[0102][0103]
式中：σ
xx
、σ
yy
、σ
zz
为x、y、z三方向主应力；τ
yz
、τ
zx
、τ
xy
为yz、zx、xy三方向切应力；ε
xx
、ε
yy
、ε
zz
为x、y、z三方向主应变；γ
yz
、γ
zx
、γ
xy
为yz、zx、xy三方向切应变；c
11-c
66
为弹性矩阵元。相应的力学参数求解公式如下(e、v公式同上)：
[0104]
体积模量b：
[0105][0106]
剪切模量g：
[0107][0108]
14)结合等效夹杂理论与有限元分析方法即可获得具有各种矿物含量的不同岩性岩石的宏观力学性质，为页岩力学性质的求取提供了一种更为简便的方法，也为目前页岩取芯困难、标准岩样制备成功率低等导致页岩力学性质测定困难等问题提供了有效解决手段。

技术特征：
1.一种基于分子动力学模拟求解矿物微观弹性参数的方法，其特征在于：从微观尺度分别建立石英、伊利石及方解石的分子动力学模型，使用msi2lmp程序将各分子动力学模型坐标转换为lammps识别的数据文件；分别描述石英内部相互作用、伊利石内部相互作用及方解石内部相互作用，通过对比驰豫前后矿物成分的晶格参数，确定所选力场合理性；通过对分子动力学模型不同方向施加形变，收集每次应变后应力应变数据，准确求解石英、伊利石及方解石的弹性常数、体积模量、剪切模量；采用delta方式对分子动力学模型模型施加应变，定义形变量为1.0
×
10-6
与施加载荷方向晶格长度之积，原子随机抖动量为分别对x、y、z、xy、xz、yz六个方向施加拉伸载荷和压缩载荷，每次操作后使用能量最小化的方式松弛原子位置以确保分子动力学模型结构稳定，其中能量的停止容差为1.0
×
10-8
，力的停止容差为1.0
×
10-10
，最大迭代次数为100，计算力或能量的最大次数为1000；不同方向进行拉伸或压缩操作后获得新的应力张量，共进行12次计算求解全部弹性常数，进一步计算石英、伊利石及方解石的体积模量、剪切模量。2.根据权利要求1所述的基于分子动力学模拟求解矿物微观弹性参数的方法，其特征在于包括如下步骤：1)分别建立石英分子动力学模型、伊利石分子动力学模型、方解石分子动力学模型，利用msilmp程序将各模型坐标转换为lammps识别的数据文件；2)描述石英原子间相互作用、伊利石原子间相互作用、方解石原子间相互作用，设置非键势能、键伸缩能、键角弯曲能以及不同原子电荷参数，各公式如下：总相互作用能：e＝e
nobonds
+e
bonds
+e
angles
非键结势能(包括l-j势能和静电势能)：键伸缩势能：键角弯曲势能：式中：ε
ij
、σ
ij
为原子i和j之间的能量参数和距离参数；r
ij
为分离距离；k
r
和k
θ
为力常数；为键的夹角；e为电子电荷；q
i
为i原子的部分电荷；q
j
为j原子的部分电荷；∈0为真空介电常数8.85419
×
10-12
f/m；r和分别为键长和键角；r0和为相应量的平衡值；ε
i
和σ
i
为常见的lennard-jones能量和尺寸参数；3)设置初始模拟参数，计量单位选择real格式，截断半径为边界条件为三维周期性边界，计算精度设为1
×
10-4
，时间步长为1fs；4)使用共轭梯度分别对石英分子动力学模型、伊利石分子动力学模型、方解石分子动力学模型进行能量最小化，对各模型初始构型进行调整；
5)定义热力学输出信息，包括步数、温度、压力、总势能、动能，晶格长度、晶格角度，设置每间隔100步输出一次；6)采用nose hoover方法控温，使各模型先后在nvt、npt系综下充分驰豫，驰豫时间为50ps，设置温度为348.73k，压力为45.82mpa，在npt系综下输出轨迹文件，每100步输出一次；7)对比驰豫前后石英、伊利石及方解石各自的晶格常数及矿物结构稳定性，验证所选力场合理；8)对石英分子动力学模型、伊利石分子动力学模型及方解石分子动力学模型执行扩胞操作，将扩胞至合适尺寸的各模型转换为数据文件；9)定义形变量为施加应变方向晶格长度
×
10-6
，原子随机抖动量为选取验证合理后的相关势函数；10)使用能量最小化的方式松弛原子位置以确保各模型结构稳定，其中能量的停止容差为1.0
×
10-8
，力的停止公差为1.0
×
10-10
，最大迭代次数为100，计算力或能量的最大次数为1000，计算各模型初始长度及变化长度，以此来定义应变；11)使用delta方式分别对各模型x、y、z、xy、xz、yz六个方向施加拉伸载荷和压缩载荷，每次仅取施加载荷方向的应变，共进行12次数据收集；12)将拉伸过程得到的弹性常数与压缩过程得到的弹性常数结果做平均处理，得到弹性常数；13)根据上述求得的弹性常数进一步求解石英、方解石和伊利石的体积模量、剪切模量，石英和方解石为三角晶系，伊利石为单斜晶系，三角晶系和单斜晶系对应的弹性矩阵和力学参数计算公式如下:三角晶系具有6个独立弹性常数，弹性范围内，其应力σ与应变ε关系通过广义胡克定律表示：相应的力学参数求解公式如下：体积模量b：b＝(2c
11
+c
33
+2c
12
+4c
13
)/9剪切模量g：g＝(7c
11
+2c
33
+12c
44-5c
12-4c
13
)/30杨氏模量e：泊松比v：
单斜晶系具有6个独立弹性常数，其应力应变关系表示为：式中：σ
xx
、σ
yy
、σ
zz
为x、y、z三方向主应力；τ
yz
、τ
zx
、τ
xy
为yz、zx、xy三方向切应力；ε
xx
、ε
yy
、ε
zz
为x、y、z三方向主应变；γ
yz
、γ
zx
、γ
xy
为yz、zx、xy三方向切应变；c
11-c
66
为弹性矩阵元；相应的力学参数求解公式如下：体积模量b：剪切模量g：单斜晶系中杨氏模量与泊松比的计算公式与三角晶系中杨氏模量与泊松比的计算公式相同。

技术总结
本发明涉及的是一种基于分子动力学模拟求解矿物微观弹性参数的方法，它包括：从微观尺度分别建立石英、伊利石及方解石的分子动力学模型，使用msi2lmp程序将各分子动力学模型坐标转换为Lammps识别的数据文件；分别描述石英、伊利石及方解石内部相互作用，通过对比驰豫前后矿物成分的晶格参数，确定所选力场合理性；通过对分子动力学模型不同方向施加形变，收集每次应变后应力应变数据，准确求解石英、伊利石及方解石的弹性常数、体积模量、剪切模量。本发明使用分子动力学模拟方法研究页岩力学性质，为目前页岩取芯困难、标准岩样制备成功率低等导致页岩力学性质测定困难，提供了有效解决手段。效解决手段。效解决手段。


技术研发人员：梁爽 孙硕 高明玉 王京坤 秦泽昊 杨公臣
受保护的技术使用者：东北石油大学
技术研发日：2023.03.05
技术公布日：2023/10/8