CN111221039B - 天然气水合物的弹性波波速确定方法、装置和设备 - Google Patents
天然气水合物的弹性波波速确定方法、装置和设备 Download PDFInfo
- Publication number
- CN111221039B CN111221039B CN202010070119.1A CN202010070119A CN111221039B CN 111221039 B CN111221039 B CN 111221039B CN 202010070119 A CN202010070119 A CN 202010070119A CN 111221039 B CN111221039 B CN 111221039B
- Authority
- CN
- China
- Prior art keywords
- crystal structure
- natural gas
- gas hydrate
- strain
- elastic
- 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.)
- Active
Links
- NMJORVOYSJLJGU-UHFFFAOYSA-N methane clathrate Chemical compound C.C.C.C.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O NMJORVOYSJLJGU-UHFFFAOYSA-N 0.000 title claims abstract description 152
- 238000000034 method Methods 0.000 title claims abstract description 60
- 239000013078 crystal Substances 0.000 claims abstract description 246
- 230000005489 elastic deformation Effects 0.000 claims abstract description 43
- 230000001902 propagating effect Effects 0.000 claims abstract description 12
- 238000000329 molecular dynamics simulation Methods 0.000 claims description 35
- 238000012545 processing Methods 0.000 claims description 30
- 238000004088 simulation Methods 0.000 claims description 29
- 230000003993 interaction Effects 0.000 claims description 12
- 230000000737 periodic effect Effects 0.000 claims description 10
- 230000009878 intermolecular interaction Effects 0.000 claims description 7
- 238000002945 steepest descent method Methods 0.000 claims description 7
- 238000011161 development Methods 0.000 abstract description 11
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Natural products C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 31
- 239000011159 matrix material Substances 0.000 description 21
- 230000008569 process Effects 0.000 description 21
- 238000010586 diagram Methods 0.000 description 11
- 239000003345 natural gas Substances 0.000 description 11
- 241001421775 Thereus Species 0.000 description 8
- -1 natural gas hydrates Chemical class 0.000 description 8
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 8
- 230000009471 action Effects 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 5
- 239000011148 porous material Substances 0.000 description 5
- 230000004044 response Effects 0.000 description 5
- 238000004590 computer program Methods 0.000 description 4
- 230000008878 coupling Effects 0.000 description 4
- 238000010168 coupling process Methods 0.000 description 4
- 238000005859 coupling reaction Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000011067 equilibration Methods 0.000 description 4
- 230000006870 function Effects 0.000 description 4
- 239000007789 gas Substances 0.000 description 4
- 238000005259 measurement Methods 0.000 description 4
- CURLTUGMZLYLDI-UHFFFAOYSA-N Carbon dioxide Chemical compound O=C=O CURLTUGMZLYLDI-UHFFFAOYSA-N 0.000 description 3
- 238000004891 communication Methods 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 229910002092 carbon dioxide Inorganic materials 0.000 description 2
- 230000003247 decreasing effect Effects 0.000 description 2
- 230000005684 electric field Effects 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 239000012535 impurity Substances 0.000 description 2
- 150000002500 ions Chemical class 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 239000001569 carbon dioxide Substances 0.000 description 1
- 230000001364 causal effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000007795 chemical reaction product Substances 0.000 description 1
- 238000002485 combustion reaction Methods 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000002803 fossil fuel Substances 0.000 description 1
- 150000004677 hydrates Chemical class 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 238000005381 potential energy Methods 0.000 description 1
- 239000000047 product Substances 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000003362 replicative effect Effects 0.000 description 1
- 239000013049 sediment Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本申请提供了一种天然气水合物的弹性波波速确定方法、装置和设备，其中，该方法包括：获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构；在目标温度和目标压强下，对第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和密度；在目标温度下、在弹性变形区间内分别向第二晶体结构施加轴向应变和剪切应变，得到第二晶体结构的应力应变关系；根据应力应变关系，确定弹性常数；根据弹性常数和第二晶体结构的密度，计算弹性波在所述第二晶体结构中传播时的波速。在本申请实施例中，可以从微观尺度准确、高效、便捷地计算天然气水合物在目标温度和目标压强下的弹性波波速，为天然气水合物的勘探和开发提供有力的理论指导。
Description
技术领域
本申请涉及水合物资源勘探与开发监测技术领域，特别涉及一种天然气水合物的弹性波波速确定方法、装置和设备。
背景技术
天然气水合物是在高压和低温的条件下，由水分子构成的笼状多面体结构和被束缚在该结构中适当尺寸的气体分子所组成的晶体化合物，天然气水合物是十分重要的“清洁能源”。由于地层孔隙介质中有天然气水合物存在的情况下，会增加孔隙介质整体的刚性和硬度，使得弹性波在该孔隙介质中传播时候，波速会明显增大，因此，在勘探或者开发自然环境中赋存的天然气水合物时候，主要是通过含天然气水合物的孔隙介质的地球物理响应，来确定天然气水合物存在与否。天然气水合物中弹性波的传播速度是非常重要的物理属性，在天然气水合物资源的勘探和开发过程中扮演着十分重要的角色，并且对天然气水合物的弹性波波速的准确识别直接影响着勘探和开发过程中的工程实践效率。
现有技术中，通常是在实验室条件下测量天然气水合物弹性波波速，具体的，是在人工合成的天然气水合物样品的一端放置弹性波发生装置，另一端放置弹性波接受器，测量弹性波从发生装置开始通过天然气水合物样品进行传播后，到达另一侧弹性波接受装置的传播时间，在已知天然气水合物样品长度的情况下，可以计算出其弹性波波速。由于在实验室合成的天然气水合物中经常伴随着残余水或者剩余气体的情况(掺有杂质)，并且形成的水合物不是一块完整的晶体结构(有晶界存在，是多晶状态)，在一些情况下还会存在多种类型天然气水合物共存的情况，即实验室合成的天然气水合物无法准确地描述实际的天然气水合物。
进一步的，由于天然气水合物样品在弹性波波速的测量过程中需要是恒温和恒压的，而保持温度和压强的过程需要人为操作来辅助完成，因此，测量过程容易受到人为因素的影响从而对测量结果产生影响。因此，采用现有技术中的技术方案无法准确、便捷地确定天然气水合物的弹性波波速。
针对上述问题，目前尚未提出有效的解决方案。
发明内容
本申请实施例提供了一种天然气水合物的弹性波波速确定方法、装置和设备，以解决现有技术中无法准确、便捷地确定天然气水合物的弹性波波速的问题。
本申请实施例提供了一种天然气水合物的弹性波波速确定方法，包括：获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构；在所述目标温度和所述目标压强下，对所述第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和所述第二晶体结构的密度；在所述目标温度下、在所述弹性变形区间内分别向所述第二晶体结构施加轴向应变和剪切应变，进行非平衡分子动力学模拟，得到所述第二晶体结构的应力应变关系；根据所述应力应变关系，确定弹性常数；根据所述弹性常数和所述第二晶体结构的密度，计算弹性波在所述第二晶体结构中传播时的波速。
在一个实施例中，在获取所述天然气水合物的第一晶体结构之前，还包括：根据天然气水合物晶胞，确定天然气水合物中各个原子的位置信息；根据所述天然气水合物中各个原子的位置信息，分别在横向、纵向和竖向上进行等量扩充，得到天然气水合物的初始晶体结构；利用最速下降法对所述初始晶体结构进行校正，得到所述天然气水合物的第一晶体结构。
在一个实施例中，在所述目标温度和所述目标压强下，对所述第一晶体结构进行平衡处理，包括：设置特征参数；控制所述第一晶体结构中所有分子的键长和键角并采用周期性边界条件控制所述第一晶体结构在横向、纵向和竖向上的边界，基于所述特征参数对所述第一晶体结构进行平衡处理。
在一个实施例中，在所述目标温度下，在所述弹性变形区间内分别向所述第二晶体结构施加轴向应变和剪切应变进行非平衡分子动力学模拟，包括：设置特征参数；基于所述特征参数，在所述弹性变形区间内向所述第二晶体结构施加轴向应变进行非平衡分子动力学模拟，其中，所述轴向应变作用于所述第二晶体结构中与横向垂直的面元上，并且所述轴向应变的方向平行于所述横向；基于所述特征参数，在所述弹性变形区间内向所述第二晶体结构施加剪切应变进行非平衡分子动力学模拟，其中，所述剪切应变作用于所述第二晶体结构中与竖向垂直的面元上，并且所述剪切应变的方向垂直于所述竖向，平行于所述横向。
在一个实施例中，所述特征参数包括以下至少之一：力场参数、目标温度、目标压强、分子之间相互作用参数、长程作用力参数、模拟时间步长和模拟时间总长度。
在一个实施例中，按照以下公式，根据所述弹性常数和所述第二晶体结构的密度，计算弹性波在所述第二晶体结构中传播时的波速：
其中，VP为弹性波在所述第二晶体结构中传播时的纵波速度；VS为弹性波在所述第二晶体结构中传播时的横波速度；K、G分别为所述第二晶体结构的体积模量和剪切模量；ρ为所述第二晶体结构的密度；C11、C12、C44为所述弹性常数；S11、S12、S44为逆弹性常数。
本申请实施例还提供了一种天然气水合物的弹性波波速确定装置，包括：获取模块，用于获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构；平衡处理模块，用于在所述目标温度和所述目标压强下，对所述第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和所述第二晶体结构的密度；处理模块，用于在所述目标温度下，在所述弹性变形区间内分别向所述第二晶体结构施加轴向应变和剪切应变进行非平衡分子动力学模拟，得到所述第二晶体结构的应力应变关系；确定模块，用于根据所述应力应变关系，确定弹性常数；计算模块，用于根据所述弹性常数和所述第二晶体结构的密度，计算弹性波在所述第二晶体结构中传播时的波速。
在一个实施例中，还包括：确定单元，用于根据天然气水合物晶胞，确定天然气水合物中各个原子的位置信息；扩充单元，用于根据所述天然气水合物中各个原子的位置信息，分别在横向、纵向和竖向上进行等量扩充，得到天然气水合物的初始晶体结构；校正单元，用于利用最速下降法对所述初始晶体结构进行校正，得到所述天然气水合物的第一晶体结构。
在一个实施例中，所述平衡处理模块包括：第一设置单元，用于设置特征参数；平衡处理单元，用于控制所述第一晶体结构中所有分子的键长和键角并采用周期性边界条件控制所述第一晶体结构在横向、纵向和竖向上的边界，基于所述特征参数对所述第一晶体结构进行平衡处理。
在一个实施例中，所述处理模块包括：第二设置单元，用于设置特征参数；第一模拟单元，用于基于所述特征参数，在所述弹性变形区间内向所述第二晶体结构施加轴向应变进行非平衡分子动力学模拟，其中，所述轴向应变作用于所述第二晶体结构中与横向垂直的面元上并且所述轴向应变的方向平行于所述横向；第二模拟单元，用于基于所述特征参数，在所述弹性变形区间内向所述第二晶体结构施加剪切应变进行非平衡分子动力学模拟，其中，所述剪切应变作用于所述第二晶体结构中与竖向垂直的面元上并且所述剪切应变的方向垂直于所述竖向，平行于所述横向。
本申请实施例还提供了一种天然气水合物的弹性波波速确定设备，包括处理器以及用于存储处理器可执行指令的存储器，所述处理器执行所述指令时实现所述天然气水合物的弹性波波速确定方法的步骤。
本申请实施例还提供了一种计算机可读存储介质，其上存储有计算机指令，所述指令被执行时实现所述天然气水合物的弹性波波速确定方法的步骤。
本申请实施例提供了一种天然气水合物的弹性波波速确定方法，可以通过获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构，并在目标温度和目标压强下，对第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和第二晶体结构的密度。进一步的，可以在目标温度下、在弹性变形区间内分别向第二晶体结构施加轴向应变和剪切应变，进行非平衡分子动力学模拟，得到第二晶体结构的应力应变关系，充分考虑了形成天然气水合物所需温度和压强的影响。可以根据上述应力应变关系，确定弹性常数，并根据弹性常数和第二晶体结构的密度，计算弹性波在第二晶体结构中传播时的波速。从而可以从微观尺度计算天然气水合物的弹性波波速，并且可以准确、高效、便捷地计算天然气水合物在目标温度和目标压强下的弹性波波速，为天然气水合物的勘探和开发中的监测提供有力的理论指导，同时节省了大量的人力、财力和物力。
附图说明
此处所说明的附图用来提供对本申请的进一步理解，构成本申请的一部分，并不构成对本申请的限定。在附图中：
图1是根据本申请实施例提供的天然气水合物的弹性波波速确定方法的步骤示意图；
图2是根据本申请具体实施例提供的非平衡分子动力学模拟轴向应变的模拟结果应力σ1与应变ε1的关系示意图；
图3是根据本申请具体实施例提供的非平衡分子动力学模拟轴向应变的模拟结果应力σ2与应变ε1的关系示意图；
图4是根据本申请具体实施例提供的XX方向、ZX方向的示意图；
图5是根据本申请具体实施例提供的非平衡分子动力学模拟剪切应变的模拟结果应力σ4与应变ε4的关系示意图；
图6是根据本申请具体实施例提供的纵波速度与温度的变化关系的示意图；
图7是根据本申请具体实施例提供的横波速度与温度的变化关系的示意图；
图8是根据本申请实施例提供的天然气水合物的弹性波波速确定装置的结构示意图；
图9是根据本申请实施例提供的天然气水合物的弹性波波速确定设备的结构示意图。
具体实施方式
下面将参考若干示例性实施方式来描述本申请的原理和精神。应当理解，给出这些实施方式仅仅是为了使本领域技术人员能够更好地理解进而实现本申请，而并非以任何方式限制本申请的范围。相反，提供这些实施方式是为了使本申请公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。
本领域的技术人员知道，本申请的实施方式可以实现为一种系统、装置设备、方法或计算机程序产品。因此，本申请公开可以具体实现为以下形式，即：完全的硬件、完全的软件(包括固件、驻留软件、微代码等)，或者硬件和软件结合的形式。
虽然下文描述流程包括以特定顺序出现的多个操作，但是应该清楚了解，这些过程可以包括更多或更少的操作，这些操作可以顺序执行或并行执行(例如使用并行处理器或多线程环境)。
自然界中形成的天然气水合物中最常见的气体分子为甲烷(CH4)，故其又称CH4水合物，天然气水合物被发现广泛分布于地球上的海洋近海沉积物和高原冻土带中。根据天然气水合物晶体单元内部所含多面体种类和数量的不同，天然气水合物可以分为Ⅰ型、Ⅱ型和H型三种结构。因为天然气(CH4)是一种重要的战略能源，1体积的CH4水合物分解后，可以释放出大约164体积的CH4分子。CH4经过燃烧后的产物是水和二氧化碳(CO2)，后者又能进入水合物的笼状结构中，以CO2水合物的形式被永久封存起来，因此，天然气水合物是十分重要的“清洁能源”。
在勘探或者开发自然环境中赋存的天然气水合物时，主要是通过含天然气水合物孔隙介质的地球物理响应，来勘探或者监测天然气水合物的存在与否。其中，天然气水合物中弹性波的传播速度是地球物理响应中一个非常重要的物理属性，在天然气水合物资源的勘探和开发过程中扮演着十分重要的角色。而现有技术中通常采用在实验室条件下测量天然气水合物弹性波波速的方式来确定天然气水合物的理论弹性波波速，具体的，是在人工合成的天然气水合物样品的一端放置弹性波发生装置，另一端放置弹性波接受器，测量弹性波从发生装置开始通过天然气水合物样品进行传播后，到达另一侧弹性波接受装置的传播时间，在已知天然气水合物样品长度的情况下，可以计算出其弹性波波速。
而采用上述方式测量得到的弹性波波速可能会受到很多因素的影响，从而使得测量结果不准确。首先，在实验室合成的天然气水合物中经常伴随着残余水或者剩余气体的情况(掺有杂质)，并且形成的水合物不是一块完整的晶体结构(有晶界存在，是多晶状态)，在一些情况下还会存在多种类型天然气水合物(例如Ⅰ型和Ⅱ型天然水合物)共存的情况，即实验室合成的天然气水合物无法准确地描述实际的天然气水合物。其次，由于天然气水合物样品的弹性波波速在测量过程中需要保证是恒温、恒压的，而保持温度和压强的过程需要人为操作来辅助完成，因此，测量过程容易受到人为因素的影响从而对测量结果产生影响。进一步的，在实验室合成的天然气水合物需要耗费一定的时间和材料，因此，采用上述方式无法高效、准确、便捷地确定天然气水合物的弹性波波速。
基于以上问题，本发明实施例提供了一种天然气水合物的弹性波波速确定方法，如图1所示，可以包括以下步骤：
S101：获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构。
由于天然气水合物在不同温度、压强条件下的弹性波波速会存在一定的差异，因此，可以预先设定需要测量在什么温度和压强条件下的天然气水合物的弹性波波速。在一个实施例中，可以预先获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构。其中，上述目标温度和目标压强可以是预先设定好的，上述目标温度可以为233K(开尔文)、238K、243K等，当然还可以采用摄氏温度来表示，温度的具体数值可以根据实际情况确定，本申请对此不作限定。
上述天然气水合物的弹性变形区间可以用于表示天然气水合物受外力作用而使各点间相对位置改变(物体发生形变)，当外力撤销后水合物可恢复原状的变形区间，对应的当外力去除时不能恢复的形变称为塑性变形。上述弹性变形区间也可称为弹性限度，是指物体受到外力作用，在内部所产生的抵抗外力的作用力不超过某一极限值时，若外力作用停止，其形变可全部消失而恢复原状。
晶体结构可以用于表示晶体的微观结构，是指晶体中实际质点(原子、离子或分子)的具体排列情况，可以根据晶体内部原子、离子、分子在空间内的三维周期性的排列规则确定。上述第一晶体结构可以为一个立方体结构，即在横向、纵向、竖向三个方向上排列的天然气水合物晶胞的数量是相同的。当然上述第一晶体结构还可以为其它形式，例如长方体、球形等，具体的可以根据实际情况确定，本申请对此不作限定。
上述天然气水合物的第一晶体结构可以为需要研究的天然气水合物的晶体结构，例如Ⅰ型天然气水合物的晶体结构、Ⅱ型天然气水合物的晶体结构或者H型天然气水合物的晶体结构等，具体的可以根据实际情况确定，本申请对此不作限定。由于天然气是一种主要由甲烷组成的气态化石燃料。因此，在一个实施例中，上述天然气水合物的晶体结构也可以是甲烷水合物的晶体结构。
S102：在目标温度和目标压强下，对第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和第二晶体结构的密度。
由于自然环境中存在的天然气水合物的晶体结构均是平衡稳定的，而上述第一晶体结构在一些情况下可能不是稳定的，因此，在一个实施例中，可以在目标温度和目标压强下，对上述第一晶体进行平衡处理，从而得到平衡后的第二晶体结构和第二晶体结构的密度。
平衡处理是指在不施加外力做功情况下，在目标温度和目标压强条件下使天然气水合物的第一晶体结构从之前非平衡状态达到该温压条件下的平衡状态，是一个能量最小化的过程。在一些实施例中，可以利用牛顿第二定律进行平衡处理，当然可以理解的是还可以采用其它任何可能的方式进行平衡处理，具体的可以根据实际情况确定，本申请对此不作限定。
在经过平衡处理之后，天然气水合物的晶体结构中含有的所有原子坐标位置都会更新，从而使整个体系的总势能降到最低。上述平衡后的天然气水合物晶体结构中可以包括以下至少之一：所有原子的坐标、分子瞬时速度、模拟体系的尺寸大小。
S103：在目标温度下、在弹性变形区间内分别向第二晶体结构施加轴向应变和剪切应变，进行非平衡分子动力学模拟，得到第二晶体结构的应力应变关系。
在一个实施例中，可以通过向第二晶体结构施加恒定速率的应变，进行非平衡分子动力学模拟(NEMD，Non equilibrium molecular dynamics)，以得到第二晶体结构在目标温度下、在弹性变形区间内产生形变的过程。从而可以得到在该形变过程中，第二晶体结构内部所产生的应力-应变关系。其中，上述恒定速率的应变可以包括：轴向应变和剪切应变。在弹性区间内，应力和应变二者之间呈线性关系。
上述应力是指当构件受到外力(也可以是温度场、电场等外因)作用时，构件内各部分之间产生相互作用的内力，用以抵抗这种外因的作用，并试图使物体从变形后的位置恢复到变形前的位置。在所考察的某个截面上，单位面积上内力称为应力，应力可以分为：拉应力、压应力、剪切应力。
上述应变是指当构件受到外力(也可以是温度场、电场等外因)作用时，几何形状和尺寸发生相对改变的物理量。应变是一个比值，没有单位。如果在应变前面加上负号，表示尺寸相对减小了，对应的，应变可以分为：拉应变、压应变、剪切应变。其中，轴向应变是指应变方向与坐标轴(X，Y或Z)相平行的应变，或者与施加的载荷在同一轴上的应变。剪切应变是剪切时物体所产生的相对形变量，即指在简单剪切的情况下，受到的力是与截面相平行的大小相等、方向相反的两个力。
由于上述第二晶体结构可以包括多个面元，每个面元上都可以施加应力从而产生应变，因此，在一个实施例中，上述轴向应变可以作用于第二晶体结构中与横向垂直的面元上，并且轴向应变的方向平行于横向，上述剪切应变可以作用于第二晶体结构中与竖向垂直的面元上，并且剪切应变的方向垂直于竖向，平行于所述横向。可以理解的是，上述轴向应变和剪切应变还可以作用于其它面元上，例如：与竖向垂直的面元上，具体的可以根据实际情况确定，本申请对此不作限定。
S104：根据应力应变关系，确定弹性常数。
由于在弹性变形区间内应力与应变之间呈线性关系，因此，可以根据上述第二晶体结构的应力-应变关系，确定弹性常数，上述弹性常数可以用于表征第二晶体结构弹性模量。
由于在弹性变形区间内物体的形变跟引起形变的外力成正比，因此，在一个实施例中，可以根据应力与应变的比值确定上述弹性常数。例如：可以按照以下公式计算弹性常数：
其中，上述C11为弹性常数；上述σ1应力；上述ε1为应变，下角标代表在矩阵中的位置。
在一个实施例中，在上述第二晶体为立方体结构的情况下，上述应力-应变关系可以如下式中所示：
其中，上述σ1为应力矩阵中第一行中的应力；C11为弹性矩阵中第一排第一列中的弹性参数；ε1为应变矩阵中第一行中的应变；n为第二晶体结构中每排每列中天然气水合物晶胞的数量，下角标用于标识各参数在矩阵中的位置，矩阵中其它参数的含义可以参照σ1、C11、ε1，此处不再赘述。
S105：根据弹性常数和第二晶体结构的密度，计算弹性波在第二晶体结构中传播时的波速。
在确定了上述第二晶体结构的弹性常数和密度之后，可以根据上述弹性常数和密度计算弹性波在第二晶体结构中传播时的波速，上述波速可以包括：横波速度和纵波速度。
在计算得到上述理论波速之后，可以将实际勘探开发过程中目标地层的地球物理响应中的弹性波波速与上述理论波速进行比对，如果实际勘探开发过程中目标地层的地球物理响应中的弹性波波速与上述理论波速相等或者比较接近，则说明目标地层中极有可能蕴含天然气水合物，从而为天然气水合物的勘探和开发提供了有力的证据，其中，实际勘探开发过程中目标地层的温度和压强与上述目标温度和目标压强相同。
在上述第二晶体结构为立方体结构，并且在横向、纵向、竖向三个方向上排列的天然气水合物晶胞的数量均为6的情况下，由于第二晶体结构属于立方晶系，所有弹性矩阵中只含有3个弹性常数：C11、C12、C44，上述应力-应变关系可以如下式中所示：
其中，上述σ1为应力矩阵中第一行中的应力；C11为弹性矩阵中第一排第一列中的弹性参数；ε1为应变矩阵中第一行中的应变，下角标用于标识各参数在矩阵中的位置，矩阵中其它参数的含义可以参照σ1、C11、ε1，此处不再赘述。
在一个实施例中，可以确定上述弹性矩阵C的逆矩阵S，从而可以求解得到逆弹性常数：S11、S12、S44。进一步的，可以根据上述弹性常数C11、C12、C44和逆弹性常数S11、S12、S44，计算上述第二晶体结构的弹性模量。在计算弹性模量时可以通过求解Voigt模型和Reuss模型二者的平均值来确定上述第二晶体结构的弹性模量，具体的，可以按照以下公式计算第二晶体结构的弹性模量：
在Voigt模型(简单线性粘弹性的力学模型)中：
其中，上述Kv为通过Voigt模型得到的体积模量；上述Gv为通过Voigt模型得到的剪切模量。
在Reuss模型(理想弹塑性模型)中：
其中，上述KR为通过Reuss模型得到的体积模量；上述GR为通过Reuss模型得到的剪切模量。
其中，VP为弹性波在第二晶体结构中传播时的纵波速度；VS为弹性波在第二晶体结构中传播时的横波速度；K、G分别为第二晶体结构的体积模量和剪切模量；ρ为第二晶体结构的密度；C11、C12、C44为弹性常数；S11、S12、S44为逆弹性常数。
在一些实施例中，上述第一晶体结构可以按照以下方式确定，具体的，可以根据天然气水合物晶胞，确定天然气水合物中各个原子的位置信息，其中，上述位置信息可以为坐标数据，原点坐标点可以为晶胞的几何中心点或者为最左下角的原子，具体的可以根据实际情况确定，本申请对此不作限定。
由于天然气水合物的晶体结构是通过在空间上复制其基本单元体(晶胞)而得到的，一个天然气水合物晶胞中包含两种由水分子组成的多面体结构，二者在晶胞中的数量不同。因此，可以根据天然气水合物中各个原子的位置信息，分别在横向、纵向和竖向上进行等量扩充，得到天然气水合物的初始晶体结构。在进行分子动力学模拟时所用模型如果太小(比如一个晶胞)，则计算结果不具有统计性，如果尺寸太大，则会导致计算成本过大、计算时间较长，因此需要综合考虑上述两个因此确定扩充的大小。在一个具体的实施例中，可以进行6×6×6的扩充，当然还可以采用其它大小的等量扩充，例如：7×7×7等，具体的可以根据实际情况确定，本申请对此不作限定。
在一些实施例中，由于得到的天然气水合物的初始晶体结构中可能会存在相邻分子空间排布过于紧密，导致相互间作用力(过剩能量)过大情况，因此，在对上述初始晶体结构进行平衡处理之前还需要对上述初始晶体结构进行校正，得到天然气水合物的第一晶体结构，以校正初始晶体结构中各个分子之间距离，使其处于合理范围之内。在一些实施例中，可以采用最速下降法或者梯度算法对上述初始晶体结构进行校正。
在一个实施例中，可以采用以下方式在目标温度和目标压强下，对第一晶体结构进行平衡处理，具体的，可以设置平衡处理的特征参数，上述特征参数可以用于表征进行平衡时的环境参数。进一步的，可以控制第一晶体结构中所有分子的键长和键角并采用周期性边界条件控制第一晶体结构在横向、纵向和竖向上的边界，基于特征参数对第一晶体结构进行平衡处理。其中，上述特征参数可以包括以下至少之一：力场参数、目标温度、目标压强、分子之间相互作用参数、长程作用力参数、模拟时间步长和模拟时间总长度。
在一个实施例中，可以采用以下方式在目标温度下，在弹性变形区间内分别向第二晶体结构施加轴向应变和剪切应变进行非平衡分子动力学模拟，具体的，可以设置非平衡分子动力学模拟的特征参数，其中，上述特征参数可以用于表征进行非平衡分子动力学模拟时的环境参数，上述特征参数可以包括但不限于以下至少之一：力场参数、目标温度、分子之间相互作用参数、长程作用力参数、模拟时间步长和模拟时间总长度。上述特征参数的具体数值可以根据实际情况确定，本申请对此不作限定。
进一步的，可以基于特征参数，在弹性变形区间内分别向第二晶体结构施加轴向应变和剪切应变进行非平衡分子动力学模拟。向第二晶体结构施加轴向应变进行非平衡分子动力学模拟，可以根据得到的模拟结果可以确定上述弹性常数C11、C12；向第二晶体结构施加剪切应变进行非平衡分子动力学模拟，可以根据得到的模拟结果可以确定上述弹性常数C44。其中，轴向应变作用于第二晶体结构中与横向垂直的面元上，并且轴向应变的方向平行于横向，剪切应变作用于第二晶体结构中与竖向垂直的面元上，并且剪切应变的方向垂直于竖向，平行于所述横向。
从以上的描述中，可以看出，本申请实施例实现了如下技术效果：可以通过获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构，并在目标温度和目标压强下，对第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和第二晶体结构的密度。进一步的，可以在目标温度下、在弹性变形区间内分别向第二晶体结构施加轴向应变和剪切应变，进行非平衡分子动力学模拟，得到第二晶体结构的应力应变关系，充分考虑了形成天然气水合物所需温度和压强的影响。可以根据上述应力应变关系，确定弹性常数，并根据弹性常数和第二晶体结构的密度，计算弹性波在第二晶体结构中传播时的波速。从而可以从微观尺度计算天然气水合物的弹性波波速，并且可以准确、高效、便捷地计算天然气水合物在目标温度和目标压强下的弹性波波速，为天然气水合物的勘探和开发中的监测提供有力的理论指导，同时节省了大量的人力、财力和物力。
下面结合一个具体实施例对上述方法进行说明，然而，值得注意的是，该具体实施例仅是为了更好地说明本申请，并不构成对本申请的不当限定。
本发明实施提供了一种天然气水合物的弹性波波速确定方法，可以包括以下步骤：
步骤1：建立甲烷(CH4)水合物晶胞(1×1×1)，并将得到的甲烷水合物晶胞结构的坐标信息转换为坐标文件，其中，坐标文件中包含甲烷水合物晶胞中所有原子的坐标位置。
步骤2：根据得到的甲烷水合物晶胞的坐标文件，在三维空间的X轴、Y轴和Z轴方向上分别进行6×6×6扩充，得到尺寸为7.218nm×7.218nm×7.218nm的甲烷水合物的初始晶体结构模型。
采用最速下降法对上述甲烷水合物的初始晶体结构模型进行校正，以校正初始晶体结构中各个分子之间距离，使其处于合理范围之内。
步骤3：设置特征参数，上述特征参数包括：力场参数、温度参数、压强参数、分子之间相互作用参数、长程作用力参数、模拟时间步长和模拟时间总长度。其中：甲烷分子力场采用OPLS-AA模型；水分子采用TIP4P/Ice模型；采用Nose-Hoover耦合控制温度参数，设置为233K，时间常数为0.5ps；采用Parrinello-Rahman耦合控制压强参数，设置为20MPa，时间常数为0.5ps；相同分子的相互作用采用Lennard-Jones参数，不同分子之间相互作用采用Lorentz-Berthelot混合规则计算；长程力作用采用Particle-Mesh-Ewald(PME)方法计算得到；模拟时间步长为时1fs；平衡时间为1ns。
此外，采用Lincs方法控制所有分子的键长和键角，并采用周期性边界条件控制校正后的初始晶体结构在X轴、Y轴和Z轴三个方向上的边界，进行平衡处理，得到平衡后的甲烷水合物的第二晶体结构(包括平衡后的原子坐标，分子瞬时速度，模拟体系大小等信息)。并可以获得在温度为233K、压强为20MPa的条件下，甲烷水合物的密度信息。
步骤4：设置特征参数，上述特征参数包括：力场参数、温度参数、分子之间相互作用参数、长程作用力参数、模拟时间步长、应变大小、应变速率。其中：甲烷分子力场采用OPLS-AA模型；水分子采用TIP4P/Ice模型；采用Nose-Hoover耦合控制温度参数，设置为233K，时间常数为0.5ps；相同分子的相互作用采用Lennard-Jones参数，不同分子之间相互作用采用Lorentz-Berthelot混合规则计算得到；长程力作用采用Particle-Mesh-Ewald(PME)方法计算得到；模拟时间步长为1fs；对平衡后的甲烷水合物晶体结构，在XX方向上施加轴向应变，其大小设置为0.04，应变速率设置为4×108s-1，其中，s代表时间单位“秒”，s-1的意思是每秒(即：/秒)。
进一步的，采用Lincs方法控制所有分子的键长和键角；采用周期性边界条件控制X轴、Y轴和Z轴三个方向上的边界，进行非平衡分子动力学模拟。由于在第二晶体上施加的轴向应变是一个递增的过程，即应变是从0到0.04根据应变速率逐渐递增的多个数值。应变速率是固定不变的，随着模拟时间的变化，应变量会逐渐增加(直到0.04为止)，在这个过程中每个应变量都会对应一个应力。因此，模拟结果可以如图2和图3中所示，其中，根据图2中的数据可以确定弹性常数C11，即图2中拟合得到的直线的斜率；根据图3中的数据可以确定弹性常数C12，即图3中拟合得到的直线的斜率。
步骤5：设置特征参数，上述特征参数包括：力场参数、温度参数、分子之间相互作用参数、长程作用力参数、模拟时间步长、应变大小、应变速率。其中：甲烷分子力场采用OPLS-AA模型；水分子采用TIP4P/Ice模型；采用Nose-Hoover耦合控制温度参数，设置为233K，时间常数为0.5ps；相同分子的相互作用采用Lennard-Jones参数，不同分子之间相互作用采用Lorentz-Berthelot混合规则计算得到；长程力作用采用Particle-Mesh-Ewald(PME)方法计算得到；模拟时间步长时1fs；对平衡后的甲烷水合物晶体结构，在ZX方向上施加剪切应变，其大小设置为0.1，应变速率设置为1×109s-1(每秒)。其中，上述XX方向、ZX方向的具体方向可以如图4中所示。
进一步的，采用Lincs方法控制所有分子的键长和键角，并采用周期性边界条件控制X轴、Y轴和Z轴三个方向上的水合物晶体结构的边界，以消除边界效应，进行非平衡分子动力学模拟。由于在第二晶体上施加的剪切应变是一个递增的过程，即应变是从0到0.1根据应变速率逐渐递增的多个数值。应变速率是固定不变的，随着模拟时间的变化，应变量会逐渐增加(直到0.1为止)，在这个过程中每个应变量都会对应一个应力。因此，模拟结果可以如图5中所示，其中，根据图5中的数据可以确定弹性常数C44，即图5中拟合得到的直线的斜率。
步骤6：根据上述弹性常数求解体积模量K和剪切模量G。
由于第二晶体结构为立方体结构，并且在横向、纵向、竖向三个方向上排列的天然气水合物晶胞的数量均为6，因此，弹性矩阵C中只含有3个弹性常数：C11、C12、C44，上述应力-应变关系可以如下式中所示：
其中，上述σ1为应力矩阵中第一行中的应力；C11为弹性矩阵中第一排第一列中的弹性参数；ε1为应变矩阵中第一行中的应变，下角标用于标识各参数在矩阵中的位置，矩阵中其它参数的含义可以参照σ1、C11、ε1，此处不再赘述。
进一步的，可以采用MATLAB软件，确定上述弹性矩阵C的逆矩阵S，从而求解得到弹性常数C11、C12、C44对应的逆弹性常数：S11、S12、S44。
在一个实施例中，上述步骤4、5、6可以利用GROMACS软件实现。
步骤7：计算弹性波在上述天然气水合物的第二晶体结构中传播时的纵波速度VP和横波速度VS。
首先是求Voigt模型和Reuss模型二者的平均，具体的，可以按照下式计算：
在Voigt模型(简单线性粘弹性的力学模型)中：
其中，上述Kv为通过Voigt模型得到的体积模量；上述Gv为通过Voigt模型得到的剪切模量。
在Reuss模型(理想弹塑性模型)中：
其中，上述KR为通过Reuss模型得到的体积模量；上述GR为通过Reuss模型得到的剪切模量。
在求得Voigt模型和Reuss模型二者的平均之后，可以根据第二晶体结构的弹性模量以及第二晶体结构的密度，按照下式计算弹性波在上述天然气水合物的第二晶体结构中传播时的纵波速度VP和横波速度VS：
其中，VP为弹性波在第二晶体结构中传播时的纵波速度；VS为弹性波在第二晶体结构中传播时的横波速度；K、G分别为第二晶体结构的体积模量和剪切模量；ρ为第二晶体结构的密度；C11、C12、C44为弹性常数；S11、S12、S44为逆弹性常数。
保持上述步骤中的其它参数不变，将步骤3、4、5中温度从233K到278K(间隔为5K)进行调整，分别计算得到在压强为20MPa，温度从233K到278K(间隔为5K)中每一个温度对应的纵波速度和横波速度，其中纵波速度与温度的变化关系图可以如图6中所示，横波速度与温度的变化关系图可以如图7中所示。
根据图6和图7中所示的内容，可知天然气水合物的弹性波波速在压强恒定的条件下，随着温度的升高而降低。其中，纵波速度随着温度的升高从4.13km/s降低到4.00km/s；横波速度随着温度的升高从1.83km/s降低到1.81km/s，由此可见，天然气水合物的刚性和抵抗外力的强度随着温度的升高而减小。
基于同一发明构思，本申请实施例中还提供了一种天然气水合物的弹性波波速确定装置，如下面的实施例。由于天然气水合物的弹性波波速确定装置解决问题的原理与天然气水合物的弹性波波速确定方法相似，因此天然气水合物的弹性波波速确定装置的实施可以参见天然气水合物的弹性波波速确定方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”指可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。图8是本申请实施例的天然气水合物的弹性波波速确定装置的一种结构框图，如图8所示，可以包括：获取模块901、平衡处理模块902、处理模块903、确定模块904和计算模块905，下面对该结构进行说明。
获取模块901，可以用于获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构；
平衡处理模块902，可以用于在目标温度和目标压强下，对第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和第二晶体结构的密度；
处理模块903，可以用于在目标温度下，在弹性变形区间内分别向第二晶体结构施加轴向应变和剪切应变进行非平衡分子动力学模拟，得到第二晶体结构的应力应变关系；
确定模块904，可以用于根据应力–应变关系，确定弹性常数；
计算模块905，可以用于根据弹性常数和第二晶体结构的密度，计算弹性波在第二晶体结构中传播时的波速。
在一个实施例中，上述天然气水合物的弹性波波速确定装置还可以包括：确定单元，用于根据天然气水合物晶胞，确定天然气水合物中各个原子的位置信息；扩充单元，用于根据天然气水合物中各个原子的位置信息，分别在横向、纵向和竖向上进行等量扩充，得到天然气水合物的初始晶体结构；校正单元，用于利用最速下降法对初始晶体结构进行校正，得到天然气水合物的第一晶体结构。
在一个实施例中，上述平衡处理模块902可以包括：第一设置单元，用于设置特征参数；平衡处理单元，用于控制第一晶体结构中所有分子的键长和键角并采用周期性边界条件控制第一晶体结构在横向、纵向和竖向上的边界，基于特征参数对第一晶体结构进行平衡处理。
在一个实施例中，上述处理模块903可以包括：第二设置单元，用于设置特征参数；第一模拟单元，用于基于特征参数，在弹性变形区间内向第二晶体结构施加轴向应变进行非平衡分子动力学模拟，其中，轴向应变作用于第二晶体结构中与横向垂直的面元上并且轴向应变的方向平行于横向；第二模拟单元，用于基于特征参数，在弹性变形区间内向第二晶体结构施加剪切应变进行非平衡分子动力学模拟，其中，剪切应变作用于第二晶体结构中与竖向垂直的面元上并且剪切应变的方向垂直于竖向，平行于所述横向。
本申请实施方式还提供了一种电子设备，具体可以参阅图9所示的基于本申请实施例提供的天然气水合物的弹性波波速确定方法的电子设备组成结构示意图，电子设备具体可以包括输入设备101、处理器102、存储器103。其中，输入设备101具体可以用于输入目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构。处理器102具体可以用于获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构；在目标温度和目标压强下，对第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和第二晶体结构的密度；在目标温度下、在弹性变形区间内分别向第二晶体结构施加轴向应变和剪切应变，进行非平衡分子动力学模拟，得到第二晶体结构的应力–应变关系；根据应力–应变关系，确定弹性常数；根据弹性常数和第二晶体结构的密度，计算弹性波在第二晶体结构中传播时的波速。存储器103具体可以用于存储弹性波在第二晶体结构中传播时的波速等参数。
在本实施方式中，输入设备具体可以是用户和计算机系统之间进行信息交换的主要装置之一。输入设备可以包括键盘、鼠标、摄像头、扫描仪、光笔、手写输入板、语音输入装置等；输入设备用于把原始数据和处理这些数的程序输入到计算机中。输入设备还可以获取接收其他模块、单元、设备传输过来的数据。处理器可以按任何适当的方式实现。例如，处理器可以采取例如微处理器或处理器以及存储可由该(微)处理器执行的计算机可读程序代码(例如软件或固件)的计算机可读介质、逻辑门、开关、专用集成电路(ApplicationSpecific Integrated Circuit，ASIC)、可编程逻辑控制器和嵌入微控制器的形式等等。存储器具体可以是现代信息技术中用于保存信息的记忆设备。存储器可以包括多个层次，在数字系统中，只要能保存二进制数据的都可以是存储器；在集成电路中，一个没有实物形式的具有存储功能的电路也叫存储器，如RAM、FIFO等；在系统中，具有实物形式的存储设备也叫存储器，如内存条、TF卡等。
在本实施方式中，该电子设备具体实现的功能和效果，可以与其它实施方式对照解释，在此不再赘述。
本申请实施方式中还提供了一种基于天然气水合物的弹性波波速确定方法的计算机存储介质，计算机存储介质存储有计算机程序指令，在计算机程序指令被执行时可以实现：获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构；在目标温度和目标压强下，对第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和第二晶体结构的密度；在目标温度下、在弹性变形区间内分别向第二晶体结构施加轴向应变和剪切应变，进行非平衡分子动力学模拟，得到第二晶体结构的应力应变关系；根据应力–应变关系，确定弹性常数；根据弹性常数和第二晶体结构的密度，计算弹性波在第二晶体结构中传播时的波速。
在本实施方式中，上述存储介质包括但不限于随机存取存储器(Random AccessMemory,RAM)、只读存储器(Read-Only Memory,ROM)、缓存(Cache)、硬盘(Hard DiskDrive,HDD)或者存储卡(Memory Card)。所述存储器可以用于存储计算机程序指令。网络通信单元可以是依照通信协议规定的标准设置的，用于进行网络连接通信的接口。
在本实施方式中，该计算机存储介质存储的程序指令具体实现的功能和效果，可以与其它实施方式对照解释，在此不再赘述。
显然，本领域的技术人员应该明白，上述的本申请实施例的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本申请实施例不限制于任何特定的硬件和软件结合。
虽然本申请提供了如上述实施例或流程图所述的方法操作步骤，但基于常规或者无需创造性的劳动在所述方法中可以包括更多或者更少的操作步骤。在逻辑性上不存在必要因果关系的步骤中，这些步骤的执行顺序不限于本申请实施例提供的执行顺序。所述的方法的在实际中的装置或终端产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行(例如并行处理器或者多线程处理的环境)。
应该理解，以上描述是为了进行图示说明而不是为了进行限制。通过阅读上述描述，在所提供的示例之外的许多实施方式和许多应用对本领域技术人员来说都将是显而易见的。因此，本申请的范围不应该参照上述描述来确定，而是应该参照前述权利要求以及这些权利要求所拥有的等价物的全部范围来确定。
以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请实施例可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。
Claims (10)
1.一种天然气水合物的弹性波波速确定方法，其特征在于，包括：
获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构；
在所述目标温度和所述目标压强下，对所述第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和所述第二晶体结构的密度；其中，所述第二晶体结构是在不施加外力做功的情况下，使所述第一晶体结构从非平衡状态达到在所述目标温度和目标压强条件下的平衡状态所得到的晶体结构；在所述目标温度和所述目标压强下，对所述第一晶体结构进行平衡处理，包括：设置特征参数；控制所述第一晶体结构中所有分子的键长和键角并采用周期性边界条件控制所述第一晶体结构在横向、纵向和竖向上的边界，基于所述特征参数对所述第一晶体结构进行平衡处理；
在所述目标温度下、在所述弹性变形区间内分别向所述第二晶体结构施加轴向应变和剪切应变，进行非平衡分子动力学模拟，得到所述第二晶体结构的应力应变关系；
根据所述应力应变关系，确定弹性常数；
根据所述弹性常数和所述第二晶体结构的密度，计算弹性波在所述第二晶体结构中传播时的波速。
2.根据权利要求1所述的方法，其特征在于，在获取所述天然气水合物的第一晶体结构之前，还包括：
根据天然气水合物晶胞，确定天然气水合物中各个原子的位置信息；
根据所述天然气水合物中各个原子的位置信息，分别在横向、纵向和竖向上进行等量扩充，得到天然气水合物的初始晶体结构；
利用最速下降法对所述初始晶体结构进行校正，得到所述天然气水合物的第一晶体结构。
3.根据权利要求1所述的方法，其特征在于，在所述目标温度下，在所述弹性变形区间内分别向所述第二晶体结构施加轴向应变和剪切应变进行非平衡分子动力学模拟，包括：
设置特征参数；
基于所述特征参数，在所述弹性变形区间内向所述第二晶体结构施加轴向应变进行非平衡分子动力学模拟，其中，所述轴向应变作用于所述第二晶体结构中与横向垂直的面元上，并且所述轴向应变的方向平行于所述横向；
基于所述特征参数，在所述弹性变形区间内向所述第二晶体结构施加剪切应变进行非平衡分子动力学模拟，其中，所述剪切应变作用于所述第二晶体结构中与竖向垂直的面元上，并且所述剪切应变的方向垂直于于所述竖向，平行于所述横向。
4.根据权利要求3所述的方法，其特征在于，所述特征参数包括以下至少之一：力场参数、目标温度、目标压强、分子之间相互作用参数、长程作用力参数、模拟时间步长和模拟时间总长度。
5.根据权利要求1所述的方法，其特征在于，按照以下公式，根据所述弹性常数和所述第二晶体结构的密度，计算弹性波在所述第二晶体结构中传播时的波速：
其中，VP为弹性波在所述第二晶体结构中传播时的纵波速度；VS为弹性波在所述第二晶体结构中传播时的横波速度；K、G分别为所述第二晶体结构的体积模量和剪切模量；ρ为所述第二晶体结构的密度；C11、C12、C44为所述弹性常数；S11、S12、S44为逆弹性常数。
6.一种天然气水合物的弹性波波速确定装置，其特征在于，包括：
获取模块，用于获取目标温度、目标压强、弹性变形区间和天然气水合物的第一晶体结构；
平衡处理模块，用于在所述目标温度和所述目标压强下，对所述第一晶体结构进行平衡处理，得到平衡后的第二晶体结构和所述第二晶体结构的密度；其中，所述第二晶体结构是在不施加外力做功的情况下，使所述第一晶体结构从非平衡状态达到在所述目标温度和目标压强条件下的平衡状态所得到的晶体结构；所述平衡处理模块包括：第一设置单元，用于设置特征参数；平衡处理单元，用于控制所述第一晶体结构中所有分子的键长和键角并采用周期性边界条件控制所述第一晶体结构在横向、纵向和竖向上的边界，基于所述特征参数对所述第一晶体结构进行平衡处理；
处理模块，用于在所述目标温度下，在所述弹性变形区间内分别向所述第二晶体结构施加轴向应变和剪切应变进行非平衡分子动力学模拟，得到所述第二晶体结构的应力应变关系；
确定模块，用于根据所述应力–应变关系，确定弹性常数；
计算模块，用于根据所述弹性常数和所述第二晶体结构的密度，计算弹性波在所述第二晶体结构中传播时的波速。
7.根据权利要求6所述的装置，其特征在于，还包括：
确定单元，用于根据天然气水合物晶胞，确定天然气水合物中各个原子的位置信息；
扩充单元，用于根据所述天然气水合物中各个原子的位置信息，分别在横向、纵向和竖向上进行等量扩充，得到天然气水合物的初始晶体结构；
校正单元，用于利用最速下降法对所述初始晶体结构进行校正，得到所述天然气水合物的第一晶体结构。
8.根据权利要求6所述的装置，其特征在于，所述处理模块包括：
第二设置单元，用于设置特征参数；
第一模拟单元，用于基于所述特征参数，在所述弹性变形区间内向所述第二晶体结构施加轴向应变进行非平衡分子动力学模拟，其中，所述轴向应变作用于所述第二晶体结构中与横向垂直的面元上并且所述轴向应变的方向平行于所述横向；
第二模拟单元，用于基于所述特征参数，在所述弹性变形区间内向所述第二晶体结构施加剪切应变进行非平衡分子动力学模拟，其中，所述剪切应变作用于所述第二晶体结构中与竖向垂直的面元上并且所述剪切应变的方向垂直于所述竖向，平行于所述横向。
9.一种天然气水合物的弹性波波速确定设备，其特征在于，包括处理器以及用于存储处理器可执行指令的存储器，所述处理器执行所述指令时实现权利要求1至5中任一项所述方法的步骤。
10.一种计算机可读存储介质，其特征在于，其上存储有计算机指令，所述指令被执行时实现权利要求1至5中任一项所述方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010070119.1A CN111221039B (zh) | 2020-01-21 | 2020-01-21 | 天然气水合物的弹性波波速确定方法、装置和设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010070119.1A CN111221039B (zh) | 2020-01-21 | 2020-01-21 | 天然气水合物的弹性波波速确定方法、装置和设备 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111221039A CN111221039A (zh) | 2020-06-02 |
| CN111221039B true CN111221039B (zh) | 2021-01-12 |
Family
ID=70827204
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010070119.1A Active CN111221039B (zh) | 2020-01-21 | 2020-01-21 | 天然气水合物的弹性波波速确定方法、装置和设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111221039B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112946737B (zh) * | 2021-01-20 | 2023-10-31 | 中国地质大学（北京） | 一种利用纵横波速度增量交汇图识别天然气水合物的方法 |
| CN113341464B (zh) * | 2021-06-04 | 2024-01-26 | 中国石油大学(北京) | 天然气水合物储层的识别方法、装置、设备和存储介质 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101074361A (zh) * | 2007-05-25 | 2007-11-21 | 东莞理工学院 | 高效复合型水合物抑制剂 |
| CN109725356A (zh) * | 2018-12-29 | 2019-05-07 | 中国地质调查局油气资源调查中心 | 一种天然气水合物开发模拟实验装置 |
| CN109871507A (zh) * | 2017-12-05 | 2019-06-11 | 中国矿业大学(北京) | 正交各向异性煤层裂隙绝对渗透率计算方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6994159B2 (en) * | 2003-11-04 | 2006-02-07 | Charles Wendland | System for extracting natural gas hydrate |
| CN105589111B (zh) * | 2016-02-01 | 2018-05-22 | 青岛海洋地质研究所 | 测量含水合物沉积介质地震波速与电磁衰减的装置及方法 |
-
2020
- 2020-01-21 CN CN202010070119.1A patent/CN111221039B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101074361A (zh) * | 2007-05-25 | 2007-11-21 | 东莞理工学院 | 高效复合型水合物抑制剂 |
| CN109871507A (zh) * | 2017-12-05 | 2019-06-11 | 中国矿业大学(北京) | 正交各向异性煤层裂隙绝对渗透率计算方法 |
| CN109725356A (zh) * | 2018-12-29 | 2019-05-07 | 中国地质调查局油气资源调查中心 | 一种天然气水合物开发模拟实验装置 |
Non-Patent Citations (1)
| Title |
|---|
| 天然气水合物的形成分布特征及其开发前景;吴茂炳，等;《中国石油勘探》;20030630;第8卷(第2期);75-78 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111221039A (zh) | 2020-06-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wang et al. | The modeling of crack propagation and coalescence in rocks under uniaxial compression using the novel conjugated bond-based peridynamics | |
| Singh et al. | A fracture-controlled path-following technique for phase-field modeling of brittle fracture | |
| Zhang et al. | A fictitious crack XFEM with two new solution algorithms for cohesive crack growth modeling in concrete structures | |
| Jansen et al. | Robust topology optimization of structures with imperfect geometry based on geometric nonlinear analysis | |
| He et al. | A thermodynamically consistent nonlocal damage model for concrete materials with unilateral effects | |
| Long et al. | Stochastic response analysis of the scaled boundary finite element method and application to probabilistic fracture mechanics | |
| CN111221039B (zh) | 天然气水合物的弹性波波速确定方法、装置和设备 | |
| Pereira et al. | A two-scale approach for the analysis of propagating three-dimensional fractures | |
| Scotta et al. | A scalar damage model with a shear retention factor for the analysis of reinforced concrete structures: theory and validation | |
| Yue et al. | An adaptive phase-field model based on bilinear elements for tensile-compressive-shear fracture | |
| Aubry et al. | Equilibrium shape of dislocation shear loops in anisotropic α-Fe | |
| Kaczmarczyk et al. | Energy consistent framework for continuously evolving 3D crack propagation | |
| Chockalingam et al. | Crystal plasticity with jacobian-free newton–krylov | |
| Krajcinovic et al. | Statistical damage mechanics—part I: theory | |
| Chen et al. | Seismogenic potential of a gouge‐filled fault and the criterion for its slip stability: Constraints from a microphysical model | |
| Xu et al. | Numerical implementation of a bounding surface plasticity model for sand under high strain-rate loadings in LS-DYNA | |
| Hashim et al. | An implicit non-ordinary state-based peridynamics with stabilised correspondence material model for finite deformation analysis | |
| Kumar et al. | An assessment of numerical techniques to find energy‐minimizing microstructures associated with nonconvex potentials | |
| Xue et al. | Fast methods for determining instabilities of elastic–plastic damage models through closed‐form expressions | |
| Shen et al. | Multivariate uncertainty analysis of fracture problems through model order reduction accelerated SBFEM | |
| Jin et al. | Elastic microplane formulation for transversely isotropic materials | |
| Cheng et al. | New technique for frictional contact on crack slip in the extended finite-element method framework | |
| Zamani et al. | Embedded interfaces by polytope FEM | |
| Baek et al. | Multiscale dynamic fracture analysis of composite materials using adaptive microstructure modeling | |
| Žalohar | Cosserat analysis of interactions between intersecting faults; the wedge faulting |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |