一种基于黏弹性波传播分析的短岩杆黏性系数测试方法与流程

本发明涉及一种岩石黏性系数测试方法，尤其涉及一种基于黏弹性波传播分析的短岩杆黏性系数测试方法。

背景技术：

在矿山开采、隧道开挖和水电站建设等岩体工程中，凿岩爆破产生的机械振动和炸药爆炸过程中产生的爆炸应力波，都会引起围岩介质的扰动并以波的形式向外传播。围岩中存在着各种结构面和节理等，外界扰动在节理岩体中传播会产生能量衰减和波形弥散，有些含有微缺陷和水分的岩石不仅表现出弹性性质也表现出一定的黏性特性，利用等效介质的方法可以将节理岩体中应力波传播等效为应力波在黏弹性材料中传播，因此开展黏弹性波传播规律的研究是探讨应力波在节理岩体中传播进而诱发矿山动力灾害的机理的关键一环。

目前基于等效介质法的黏性系数测试主要是利用已建立的黏弹性本构关系和动态冲击试验中测得的岩杆上两测点(p1和p2)的应变波形，按照式(1)计算应力波在岩杆上传播的衰减系数，然后根据衰减系数和本构方程中黏性系数的关系进行求解。若黏弹性本构关系取kelvin-voigt模型的本构方程，黏性系数和衰减系数的关系如式(2)所示。

依据岩杆上两测点p1和p2测得的应变波进行傅里叶变换可求得各个谐波对应的衰减系数，如式(1)所示。

式中，α为衰减系数，ω为傅里叶变换后的谐波角频率，ε2和ε1分别为p2点和p1点测得的应变，f表示对应变波进行傅里叶变换，l为p1点和p2点之间的距离，re表示取复数的实部。

利用kelvin-voigt模型建立黏弹性本构方程，然后结合一维波传播理论构建一维黏弹性波传播控制方程，黏性系数可由方程(2)求得。

上述黏性系数的测试方法存在如下不足：首先，由方程(1)求得的衰减系数是频率的函数，因此最后计算得到的黏性系数也是频率的函数并非一个定值，需要依据主频确定岩杆的黏性系数；其次，测试岩杆上两点的应变波形，需要保证岩杆的长度能够测得两个完整的加载波形，对岩杆的长度要求较高，而钻取的天然岩杆长度比较短，因此满足要求的岩杆制作比较困难。

技术实现要素：

本发明实施例所要解决的技术问题在于，提出了室内试验和数值模拟相结合的试验方法，通过试错的方法不断调整黏性系数，优化出一个最合适的黏性系数。

为解决上述技术问题，提供了一种基于黏弹性波传播分析的短岩杆黏性系数测试方法，包括以下步骤：

建立岩杆的黏弹性本构方程和一维黏弹性波传播的控制方程；

基于一维黏弹性波传播的控制方程，应用透反射理论和黏弹性损伤理论设计一维黏弹性应力波分析程序；

利用一维黏弹性应力波分析程序模拟冲击试验，输入加载波并将模拟冲击试验获得的应变波与验证波进行对比，通过调整设置的参数使模拟的应变波形与验证波形一致，进而获得岩杆的黏性系数。

在本发明基于黏弹性波传播分析的短岩杆黏性系数测试方法中，所述建立岩杆的黏弹性本构方程具体为：

利用摆锤加载中应变率霍普金森压杆试验系统进行岩杆的shpb试验，测试其动态应力-应变曲线；

借助弹性元件和黏性元件构建本构模型，结合岩杆的动态应力-应变曲线建立岩杆的动态本构方程。

在本发明基于黏弹性波传播分析的短岩杆黏性系数测试方法中，所述建立一维黏弹性波传播的控制方程具体为：

基于一维波传播的连续性方程、运动方程和岩杆的动态本构方程，建立一维黏弹性波传播的控制方程。

在本发明基于黏弹性波传播分析的短岩杆黏性系数测试方法中，利用一维黏弹性应力波分析程序模拟冲击试验的过程包括以下步骤：

前处理步骤：建立岩杆试样的几何模型，进行网格划分将模型划分成多个单元，对各单元的相关参数进行赋值，设置边界条件和加载条件；

应力分析步骤：进行节点的应力计算同时判断单元是否发生损伤并对损伤情况进行处理；

后处理步骤：输出关键节点的应力值、应变值以及弹性模量值。

在本发明基于黏弹性波传播分析的短岩杆黏性系数测试方法中，在所述前处理步骤中对以下参数进行赋值：密度、弹性模量、横截面积、黏性系数、单轴抗压强度以及单轴抗拉强度；

在前处理步骤中设置加载条件具体为：对施加在岩杆试样上的加载形式、加载波形以及波形参数进行设置。

在本发明基于黏弹性波传播分析的短岩杆黏性系数测试方法中，应力分析步骤具体包括以下步骤：

对岩杆试样施加一个新的边界应力；

根据控制方程和透反射理论计算节点的应力值和应变值；

根据破坏准则判断是否有单元损伤，若有损伤，则根据黏弹性损伤理论计算损伤情况，并修改相关物理参数后重新计算各节点的应力值和应变值；

若没有损伤，则输出各节点的应力值和应变值。

在本发明基于黏弹性波传播分析的短岩杆黏性系数测试方法中，黏性损伤根据如下公式进行计算：

η＝(1-d)η0；

其中，

式中d是损伤因子，处于0到1之间，η和η0分别是损伤后和损伤前的黏性系数；

ε为单元的应变，εt0为弹性极限应变，εtu为单元的拉伸极限应变，εc0为弹性极限对应的压应变；

e0损伤前的弹性模量，ftr为拉伸破坏残余强度，fcr为压缩破坏残余强度。

在本发明基于黏弹性波传播分析的短岩杆黏性系数测试方法中，后处理步骤具体为：

从应力分析步骤的计算结果中提取关键节点的应力值、应变值、弹性模量值，并输出应力图、应变图以及弹性模量图。

在本发明天然短岩杆黏性系数测试方法中，根据关键节点的应力值、应变值、弹性模量值绘制应变波形图，将应变波与验证波进行对比，通过调整参数使应变波形与验证波形相一致，进而获得岩杆的黏性系数。

在本发明天然短岩杆黏性系数测试方法中，利用摆锤加载中应变率霍普金森压杆试验系统进行天然短岩杆的冲击试验，获得两个测点的波形分别作为模拟冲击试验中所需的加载波形和验证波形。

上述技术方案至少具有如下有益效果：本发明提出了室内试验和数值模拟相结合的方法，通过shpb试验测得岩石的动态本构曲线进而构建相应的数值模型，然后利用动态冲击试验测得数值模拟中需要的加载波形和验证波形，通过试错的方法不断调整黏性系数，优化出一个最合适的黏性系数。

针对现有试验方法的不足，本发明的试验方法有如下改进：首先，通过试验测得的动态应力-应变曲线构建动态本构关系更加符合岩石的动态特性，测得岩石的黏性系数更加合理；其次，通过数值模拟不断试错的方法得到的黏性系数相比主频对应黏性系数更加具有代表性；再次，本发明的动态冲击试验中只需要测得一个完整的波形作为数值模型的加载波形，另一个验证波形并不要求是否有叠加和完整，因此可以缩短试验中岩杆的长度，有利于天然岩石黏性系数的测试。

附图说明

图1是kelvin-voigt模型示意图；

图2细观单元单轴受拉时的弹性损伤本构关系示意图；

图3是一维黏弹性应力波分析程序流程图。

具体实施方式

本发明提供了一种基于黏弹性波传播分析的短岩杆黏性系数测试方法，包括以下步骤：建立岩杆的黏弹性本构方程和一维黏弹性波传播的控制方程；基于一维黏弹性波传播的控制方程、透反射理论以及黏弹性损伤理论设计一维黏弹性应力波分析程序；利用一维黏弹性应力波分析程序模拟冲击试验，输入加载波并将模拟冲击试验获得的应变波与验证波进行对比，通过调整设置的参数使模拟的应变波形与验证波形一致，进而获得岩杆的黏性系数。

其中，建立岩杆的黏弹性本构方程具体为：利用摆锤加载中应变率霍普金森压杆试验系统，进行岩杆的shpb试验，测试岩杆试样的动态应力-应变曲线。借助弹性元件和黏性元件构建本构模型，弹性元件可采用弹簧，黏性元件采用黏壶。具体实施时可采用图1所示的开尔文体元件模型作为黏弹性材料的本构模型，结合shpb试验获得的岩杆的动态应力-应变曲线建立岩杆的动态本构方程，如式(3)所示。

式中，σ、ε分别是开尔文体的应力和应变，e和η分别是开尔文体中弹簧的弹性模量和牛顿体的黏性系数。

建立岩石的动态本构方程后，基于一维波传播的连续性方程、运动方程和岩石的动态本构方程，建立一维黏弹性波传播的控制方程。

由一维纵波的连续性方程可知：

由牛顿定理可得一维纵波的运动方程：

式中，v是质点速度，ρ是介质的密度，x是质点的坐标。

根据公式(3)(4)(5)可以求得一维黏弹性波的控制方程：

方程(6)的谐波解为：

u(x,t)＝u0exp(αx)exp[i(ωt-kx)](7)

式中，α和k分别为衰减因子和波数，α为负数表示谐波随传播距离指数衰减，u和ω分别表示位移和角频率。

应力波传播的相速度c＝ω/k，是频率的函数，由方程(8)可得

基于上述过程构建一维黏弹性波传播的控制方程后，应用透反射理论和黏弹性损伤理论设计一维波传播分析程序。下面具体说明一下一维波传播分析程序的设计原理：

1、应力波在黏弹性杆界面处透反射

一维谐波在无界面黏弹性杆中传播时，谐波幅值仅随传播距离按照指数函数衰减，谐波频率和相速度不发生变化。由式(10)可知，谐波的传播速度不仅受杆的材料属性影响而且也与谐波频率有关。一维应力波在弹性杆中传播遇到材料属性变化的界面时，应力波透过界面的波形幅值依据透反射关系增减，但是频率保持不变，即应力波波形图纵向伸缩而横向不变。谐波在黏弹性杆中传播时，黏弹性杆的材料属性并不影响谐波的频率，因此当谐波在黏弹性杆中传播时遇到材料属性变化的分界面时，谐波通过分界面之后的波形频率也保持不变，即谐波波形图也是纵向伸缩而横向不变。

由连续性条件可知：

fi+fr＝ft(11)

vi+vr＝vt(12)

由动量守恒得：

fidt＝ρ1a1c1dtvifi＝σia1(13)

frdt＝-ρ1a1c1dtvrfr＝σra1(14)

ftdt＝ρ2a2c2dtvtft＝σta2(15)

式中：f、v和σ分别为分界面上的力、质点速度和应力，下标i、r和t分别表示入射波、反射波和透射波，ρ1和ρ2分别为分界面左边和右边的密度，c1和c2分别为分界面左边和右边的谐波速度，dt为单位时间，a1和a2分别为分界面左边和右边的截面积。

由式(11)-(15)得：

由式(16)和(17)可知，谐波通过材质和横截面积变化的界面时，谐波的透射系数和反射系数与界面两边的材料密度、谐波相速度和横截面面积有关。由于谐波在黏弹性杆中传播的相速度受材料属性和谐波频率的影响，因此谐波的透射系数和反射系数受谐波的频率影响，不同频率的谐波透反射系数不同。由不同频率的谐波叠加形成的应力波通过材质和横截面积变化的界面时，由于不同频率谐波的透反射系数不一致，导致叠加之后形成的应力波的频率将发生变化，即应力波通过波阻抗变化的界面后，应力波波形不仅在纵向伸缩而且在横向也有伸缩，这与任意形状的应力波在一维弹性杆中传播的透反射规律不同。

2、黏弹性波加载边界处理

由于一维黏弹性波控制方程中应变和应力方程形式一样，因此应力或应变谐波在黏弹性长杆中传播规律一致。但是由于加载边界的不一致，需要进行相应的处理。

(1)应变边界

模型的加载边界为应变控制，需要由应变求解应力，由方程(3)开尔文体的本构方程可知，基于应变波传播的质点应力求解如下：

式中ε^t、σ^t、e^t和η^t表示t时刻某质点的应变、应力、弹性模量和黏性系数，δt表示时间增量。质点某时刻的应力可以由这一时刻和上一时刻的应变进行求解。

(2)应力边界

模型的加载边界为应力加载，模型中波以应变的形式进行传播，需要把应力波转换成应变波，由方程(20)可知

应力波加载的初始时刻，t＝0，ε＝0。

3、黏性损伤理论

在一维黏弹性长杆中黏弹性材料单元所受拉应力或压应力大于单轴抗拉强度或抗压强度时材料发生损伤，损伤包括两部分，一部分是弹性损伤，通过改变材料的杨氏模量来表达损伤的影响，另一部分是黏性损伤，通过改变黏性系数来表达损伤对材料的影响。

(1)弹性损伤原理

按照弹性损伤理论，单元产生损伤后，材料参数要进行弱化，弹性模量按照式(22)进行处理。

e＝(1-d)e0(22)

式中d是损伤因子，处于0到1之间；e和e0分别是损伤后和损伤前的弹性模量。

假定细观单元的弹性损伤本构关系如图2所示。当单元所受的拉应力大于材料的抗拉强度时，单元出现拉伸损伤，此时单元的应变为ε，可知εt0≤ε<εtu，εt0为弹性极限对应的应变即弹性极限应变，εtu为单元的极限拉伸应变。由于该单元尚未完全破坏仍然有一定的承载能力，而这种承载能力可以定义一个残余强度ftr来表征，损伤因子的计算如式(23)所示。

当单元所受的压应力大于材料的抗压强度时，单元出现压损伤，采用和拉伸损伤相似的方法处理压缩损伤，如式(24)所示。

式中，εc0为弹性极限对应的压应变，fcr为压缩破坏残余强度。

(2)黏性损伤原理

弹性损伤因子的计算是利用残余强度、材料损伤前的弹性模量和损伤时的应变。当材料为黏弹性材料时，单元的应力由弹性项和黏性项两部分组成，借助于弹性材料损伤因子的计算方法来计算黏弹性材料的损伤因子。拉伸损伤和压缩损伤的损伤因子分别如下式(25)和(26)所示。

参照弹性损伤中杨氏模量的计算公式给出黏性损伤中完整岩石的黏性系数和损伤后岩石黏性系数的方程式(27)。

η＝(1-d)η0(27)

式中η和η0分别是损伤后和损伤前的黏性系数。

根据上述关于应力波在一维黏弹性杆中的传播理论、透反射理论和黏弹性损伤理论的设计一维黏弹性应力波分析程序，程序流程如图3所示。再采用一维黏弹性应力波分析程序模拟冲击试验测试的应变波，将应变波形与验证波形进行对比，利用试错的方法计算天然岩石的黏性系数。具体实施过程如下：

如图3所示，利用一维波传播分析程序模拟冲击试验过程包括前处理、应力分析和后处理三部分。前处理步骤主要功能是建立岩杆试样的几何模型，进行网格划分将模型划分成多个单元，对各单元的相关参数进行赋值，设置边界条件和加载条件。应力分析步骤主要功能是进行节点应力计算同时判断单元是否发生损伤并对损伤情况进行处理。后处理步骤主要功能是输出关键节点的应力值、应变值以及弹性模量值。

前处理步骤中，首先对各单元的以下参数进行赋值：密度、弹性模量、横截面积、黏性系数、单轴抗压强度以及单轴抗拉强度。再对施加在岩杆试样上的加载形式、加载波形以及波形参数进行设置。具体实施时，先利用摆锤加载中应变率霍普金森压杆试验系统进行天然短岩杆的冲击试验，进而获得两个测点的波形。将其中一个完整的波形作为模拟冲击试验中的加载波形和另一个波形作为验证波形。

前处理步骤完成后，进行应力分析步骤。应力分析步骤具体包括以下步骤：

对岩杆试样施加一个新的边界应力；

根据控制方程和透反射理论计算各节点的应力值和应变值；

根据破坏准则判断是否有单元损伤，若有损伤根据黏性损伤理论和弹性损伤理论计算损伤情况，并修改相关物理参数后重新计算各节点的应力值和应变值；

若没有损伤，则输出各节点的应力值和应变值。

从应力分析步骤的计算结果中提取关键节点的应力值、应变值、弹性模量值，并输出应力图、应变图以及弹性模量图。

根据关键节点的应力值、应变值、弹性模量值绘制应变波形图，将应变波与验证波进行对比，通过调整参数使应变波形与验证波形相一致，进而获得岩杆的黏性系数。

本发明采用室内试验和数值模拟相结合的测试方法，首先建立岩杆的黏弹性本构方程和一维黏弹性波传播的控制方程；基于一维黏弹性波传播的控制方程、透反射理论以及黏弹性损伤理论设计一维黏弹性应力波分析程序；利用一维黏弹性应力波分析程序模拟冲击试验，输入加载波并将模拟冲击试验获得的应变波与验证波进行对比，通过调整设置的参数使模拟的应变波形与验证波形一致，进而获得岩杆的黏性系数。

本发明通过试验测得的动态应力-应变曲线构建动态本构关系更加符合岩石的动态特性，测得岩石的黏性系数更加合理；通过数值模拟不断试错的方法得到的黏性系数相比主频对应黏性系数更加具有代表性；本发明的动态冲击试验中只需要测得一个完整的波形作为数值模型的加载波形，另一个验证波形并不要求是否有叠加和完整，因此可以缩短试验中岩杆的长度，有利于天然岩石黏性系数的测试。

以上所述是本发明的具体实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。