一种适用于地震荷载下岩石动态本构模型的构建方法与流程

1.本发明属于岩体工程技术领域，更具体地，涉及一种适用于地震荷载下岩石动态本构模型的构建方法。


背景技术：

2.岩石地下工程作为国家战略和生命线工程的重要基础设施，例如，公路隧道、铁路隧道、地下洞室等，一旦面临复杂地震地质环境，其地震安全问题突出，因此，建立适用于地震荷载作用下岩石的动态本构模型，对岩石地下工程安全稳定性评价具有重要意义。
3.地震荷载由于其复杂性，往往简化为动力荷载或者循环荷载，目前岩石在动力荷载的本构模型主要是基于应变率效应和损伤效应的模型，而在循环荷载下的本构模型则是反映滞回圈和累积变形量等力学特性的模型。基于应变率效应和损伤效应，主要为经验和半经验模型(钱七虎,戚承志.岩石、岩体的动力强度与动力破坏准则[j].同济大学学报(自然科学版),2008,36(12):1599
‑
1605.)、力学模型(谢理想,赵光明,孟祥瑞.软岩及混凝土材料损伤型黏弹性动态本构模型研究[j].岩石力学与工程学报,2013,(04):857
‑
864.)和组合模型(陆晓霞,张培源.在围压冲击条件下岩石损伤粘塑性本构关系[j].重庆大学学报(自然科学版),2002,25(1):6
‑
8.)。上述模型可以较好地反映岩石材料的率效应或损伤效应，但却不能很好地模拟岩石材料在循环荷载下的滞回圈及累积变形等现象，不能很好体现岩石材料的应力路径性质。对于反映循环荷载下岩石材料滞回圈和累积变形量等力学特性的本构模型研究，主要包括内变量疲劳本构模型(王者超,赵建纲,李术才,等.循环荷载作用下花岗岩疲劳力学性质及其本构模型[j].岩石力学与工程学报,2012,31(9):1888－1900.)、二元介质本构模型(刘恩龙,张建海,何思明,等.循环荷载作用下岩石的二元介质模型[j].重庆理工大学学报(自然科学),2013,27(9):6
‑
11,16.)、损伤模型(周永强，盛谦，付晓东，冷先伦.一种适用于循环荷载下岩石损伤本构模型的构建方法，cn201911231456.8)等。上述模型很难很好地体现对于岩石材料强度和刚度在动态荷载下率效应的强化性。
[0004]
综上所述，现在构造地震荷载下岩石动态本构模型仍存在一定的缺陷，主要体现在：
[0005]
(1)实验层面上下岩石进行宏观动态力学特性的研究较多，但荷载也仅仅局限在常规的动态荷载或者循环荷载，却很难针对于真实的地震荷载；
[0006]
(2)地震荷载不但是动态荷载，也是循环荷载，目前技术方案中并没有综合考虑循环荷载下的变形特性和动态荷载下的力学特性，从而无法精确获取地震载荷下的精确力学特性；
[0007]
(3)针对地震载荷，如何将动态荷载和循环载荷结合起来反应岩石在地震荷载下的力学特征性质和变形特征，现有技术中也是没有提出的。


技术实现要素：

[0008]
针对现有技术的以上缺陷或改进需求，本发明提供了一种适用于地震荷载下岩石动态本构模型的构建方法，其综合考虑动态荷载和循环载荷下的变形特性和动态荷载下的力学特性，并结合次加载面理论、基于修正cwfs模型，采用强度或刚度动态增强因子模型来描述岩石强度或刚度与应变率的关系，并将cwfs模型与动态增强因子模型进行耦合，同时基于岩石弹性模量的损伤和应变率效应的耦合函数，从而将动态载荷和循环载荷较好的结合，同时结合偏应力
‑
轴向应变关系曲线、应力应变曲线以及岩石循环加卸载试验应力
‑
应变曲线构建和求解弹塑性矩阵，从而建立了能反应岩石在地震荷载下的力学特征性质和变形特征的动态本构模型，同时，模型参数具有物理意义并均可基于试验结果获取，具有广泛的适用性。
[0009]
为实现上述目的，本发明提出了一种适用于地震荷载下岩石动态本构模型的构建方法，包括以下步骤：
[0010]
s1基于次加载面理论和drucker
‑
prager屈服准则，构建岩石次加载面的函数以及岩石相似中心面的函数；
[0011]
s2基于试验数据绘制岩石在不同围压下的偏应力
‑
轴向应变关系曲线、不同应变率下岩石单轴压缩的应力应变曲线以及岩石循环加卸载试验应力
‑
应变曲线；
[0012]
s3将修正cwfs模型与动态增强因子模型进行耦合，以获取岩石的粘聚力的耦合计算模型，根据该耦合计算模型以及基于岩石弹性模量的损伤和应变率效应的耦合函数，并结合偏应力
‑
轴向应变关系曲线、应力应变曲线以及岩石循环加卸载试验应力
‑
应变曲线构建和求解弹塑性矩阵；
[0013]
s4根据所述弹塑性矩阵构建岩石动态本构模型，同时更新步骤s1中岩石正常屈服面与次加载面相似比的函数式，重复步骤s1至步骤s4，直至所构建的岩石动态本构模型计算的应力满足次加载面上的屈服函数。
[0014]
作为进一步优选的，步骤s1具体包括以下步骤：
[0015]
s11选择drucker
‑
prager屈服面作为次加载面理论的正常屈服面，构建正常屈服面上的屈服函数
[0016][0017]
其中，β(κ)为第二材料参数，κ为内变量，i1和j2则分别是应力σ的第一不变量和应力的偏应力第二不变量；
[0018]
s12基于次加载面理论，构建岩石次加载面的函数
[0019][0020]
其中，为次加载面上的函数，为次加载面上应力，和则分别是次加载面上应力的第一不变量和次加载面上应力的偏应力第二不变量，r为次加载面与正常屈服面的相似比，β(κ)为第二材料参数，κ为内变量，为第一材料参数；
[0021]
s13构建岩石相似中心面的函数
[0022][0023]
其中，s为相似中心，和则分别是相似中心面上相似中心s的第一不变量和相似中心面上相似中心s的偏应力第二不变量，r
s
为相似中心面与正常屈服面之比。
[0024]
作为进一步优选的，步骤s11中，
[0025]
所述第二材料参数β(κ)的表达式为：
[0026][0027]
所述第一材料参数的表达式为：
[0028][0029]
式中，所述coh为岩石的粘聚力，该粘聚力是随内变量κ和应变率变化的函数。
[0030]
作为进一步优选的，步骤s3中，所述岩石的粘聚力的耦合计算模型为：
[0031][0032]
式中，所述为岩石的内摩擦角，该内摩擦角是随内变量κ变化的函数，其具体计算模型为：
[0033][0034]
所述内变量κ的具体计算模型为：
[0035][0036][0037]
其中，c0为粘聚力初始值，c
r
为粘聚力最终值，为内摩擦角的初始值，为内摩擦角的最终值，是粘聚力开始变化时内变量的值，是内摩擦角开始变化时内变量的值，是粘聚力达到最终值时内变量的值，是内摩擦角达到最终值时内变量的值，dε
p
是
塑性应变增量，p是围压，f
c
是单轴抗压强度，a1是第三材料参数，a2是第四材料参数，tr为矩阵的迹，为静态荷载下的应变率，d为强度动态增强因子模型的第五材料参数，n是围压p的函数，其用于表示内变量k的增长速率。
[0038]
作为进一步优选的，步骤s3中，所述单轴抗压强度f
c
根据岩石单轴压缩试验求解，根据三轴分级循环加卸载试验得出不同围压p下的应力应变曲线，把峰值前和峰值后的曲线分开考虑，在峰值前阶段，首先计算不同围压p下从初始屈服到峰值间的塑性应变ε
p
，假设内变量κ在初始屈服时为0，在达到峰值时为1，然后得出不同围压p下内变量增量dκ与围压p的关系，进而求解出n的关系式；然后将内变量κ等间隔取值，将试验得到的偏应力根据对应的内变量κ的数值进行线性插值，从而得到不同围压p下相同内变量κ对应的偏应力，根据摩尔圆即可求出不同内变量κ对应的强度参数粘聚力coh和内摩擦角接着，求出峰后阶段求出不同内变量κ对应的强度参数粘聚力coh和内摩擦角根据峰值前的塑性应变ε
p
与峰后的塑性应变ε
p
的比值，把峰值前和峰后的内变量κ进行统一，从而得出强度参数粘聚力coh和内摩擦角与内变量κ的关系，进而可以求解c0，c
r
，a1，a2，根据不同应变率下岩石单轴压缩的应力应变曲线，获取不同应变率下的岩石的峰值强度，通过强度增强因子模型拟合得出参数d。
[0039]
作为进一步优选的，所述步骤s3具体包括以下步骤：
[0040]
s31求解正常屈服面上背应力α的增量dα：
[0041][0042]
式中，i为第六材料参数，r为第七材料参数，||||表示向量的模，dε
p
为应力σ对应的塑性应变增量；
[0043]
s32获取应力σ与相似中心s的矢量差次加载面上考虑正常屈服面上背应力α的应力以及相似中心s与正常屈服面上背应力α的矢量差
[0044][0045][0046][0047]
式中，σ
y
是应力σ在正常屈服面上的对偶应力，为相似中心s与次加载面上背应力的矢量差；
[0048]
s33构建并求解弹塑性矩阵d
ep
，该弹塑性矩阵d
ep
的计算模型如下：
[0049][0050]
式中，
[0051][0052][0053][0054][0055][0056][0057]
其中，“·”为点乘，d
el
为弹性矩阵，ν为泊松比，χ为相似中心比r
s
的最大值，δ是二阶单位矩阵，i是对称的四阶单位矩阵，m和u分别为第八材料参数和第九材料参数，a0和b分别为第十材料参数和第十一材料参数，是弹性模量，是随内变量κ和应变率变化的函数，e
s
和e
e
分别是粘聚力初始值和最终值，κ
e
是弹性模量达到最终值时内变量的值。
[0058]
作为进一步优选的，步骤s3中，基于试验数据绘制岩石在不同围压下的偏应力
‑
轴向应变关系曲线，根据初始阶段偏应力
‑
轴向应变关系曲线直线段，计算得出弹模初始值e
s
、泊松比ν，根据残余阶段偏应力
‑
轴向应变关系曲线直线段，计算得出弹模最终值e
e
；
[0059]
根据不同应变率下岩石单轴压缩的应力应变曲线，获取不同应变率下的岩石的弹模初始值，通过刚度增强因子模型拟合得出中的参数a0和b；
[0060]
u为控制应力点向正常屈服状态靠近的速率的参数，通过屈服状态下岩石循环加卸载试验应力
‑
应变曲线的斜率确定。
[0061]
作为进一步优选的，所述岩石正常屈服面与次加载面相似比的函数式中，将更新后的次加载面与正常屈服面的相似比r更新为r1，所述相似比r1的计算公式如下：
[0062][0063]
其中，
[0064][0065][0066][0067]
式中，为的平均值，为的平均值，为的偏应力，为的偏应力，β(κ)为第二材料参数，κ为内变量，为第一材料参数，为应力σ与相似中心s的矢量差，为相似中心s与正常屈服面上背应力α的矢量差。
[0068]
作为进一步优选的，所述岩石动态本构模型的计算公式如下：
[0069][0070]
其中，dσ为应力增量，dε为应变增量，为次加载面上的函数。
[0071]
作为进一步优选的，基于有限元软件，对不同加载速率循环荷载试验进行数值模拟，验证所述岩石动态本构模型的适用性和准确性。
[0072]
总体而言，通过本发明所构思的以上技术方案与现有技术相比，主要具备以下的技术优点：
[0073]
1.本发明综合考虑动态荷载和循环载荷下的变形特性和动态荷载下的力学特性，基于次加载面理论和修正cwfs模型，采用强度或刚度动态增强因子模型来描述岩石强度或刚度与应变率的关系，并将cwfs模型与动态增强因子模型进行耦合，同时基于岩石弹性模量的损伤和应变率效应的耦合函数以及均可基于试验结果获取的具有物理意义模型参数一同构建能反应岩石在地震荷载下的力学特征性质和变形特征的动态本构模型，使得本发明提供的动态本构模型能更真实的反应地震载荷下岩石的力学行为，使得岩石损伤模拟更加精准。
[0074]
2.本发明克服传统模型只能模拟常规的动态荷载或者循环荷载的问题，通过采用强度或刚度动态增强因子模型来描述岩石强度或刚度与应变率的关系，并将cwfs模型与动态增强因子模型进行耦合，同时，将岩石弹性模量的损伤和应变率效应进行耦合，从而将动态载荷和循环载荷较好的结合，以综合考虑循环荷载下的变形特性和动态荷载下的力学特性，更精确和真实的反应地震载荷下岩石的力学行为。
[0075]
3.本发明建立的动态本构模型参数物理意思明确，均开通通过实验试验结果获
取，且可变异为有限元等嵌入，方法简单便捷，适用性强，易于推广应用于综合考虑动态荷载和循环载荷下的变形特性和动态荷载下的力学特性的实际岩石工程计算与分析。
附图说明
[0076]
图1是本发明优选实施例涉及的一种适用于地震荷载下岩石动态本构模型的构建方法的流程图；
[0077]
图2是t
2y
6大理岩在常规三轴压缩下轴向偏应力
‑
轴向应变关系曲线、轴向偏应力
‑
径向应变关系曲线；
[0078]
图3是t
2y
6大理岩在围压为0mpa的循环荷载下轴向偏应力
‑
轴向应变关系曲线、轴向偏应力
‑
径向应变关系曲线；
[0079]
图4是t
2y
6大理岩在围压为10mpa的循环荷载下轴向偏应力
‑
轴向应变关系曲线、轴向偏应力
‑
径向应变关系曲线；
[0080]
图5是t
2y
6大理岩在围压为40mpa的循环荷载下轴向偏应力
‑
轴向应变关系曲线、轴向偏应力
‑
径向应变关系曲线；
[0081]
图6是是t
2y
6大理岩在不同应变率循环荷载下轴向偏应力
‑
轴向应变关系曲线、轴向偏应力
‑
径向应变关系曲线(包络线)和植入地震荷载下岩石的动态本构模型的有限元商业软件获得的轴向偏应力
‑
轴向应变关系曲线、轴向偏应力
‑
径向应变关系曲线；
[0082]
图7是t
2y
6大理岩的粘聚力和内摩擦角与内变量的关系图；
[0083]
图8是t
2y
6大理岩峰值强度与应变率的关系图；
[0084]
图9是t
2y
6大理岩弹性模量与内变量的关系图；
[0085]
图10是t
2y
6大理岩弹性模量与应变率的关系图；
[0086]
图11是t
2y
6大理岩在不同应变率循环荷载下累积轴向塑性应变与相对循环次数试验曲线和植入地震荷载下岩石的动态本构模型的有限元商业软件获得的数值理论累积轴向塑性应变与相对循环次数试验曲线。
具体实施方式
[0087]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0088]
本发明提供的一种适用于地震荷载下岩石动态本构模型的构建方法的原理包括：次加载面理论是由日本学者hashiguchi提出的，其基本思路是假设在正常屈服面(常规模型的屈服面)的内部存在一个与之保持几何相似的次加载面，当前应力点始终位于该加载面上，因此加载准则不需要判断应力点是否位于屈服面上。在循环荷载下，岩石的强度会劣化，而劣化的表现主要是因为岩石的粘聚力降低，内摩擦角升高，即cwfs模型。本发明通过岩石材料循环荷载试验，以不同围压下等效塑性剪应变为内变量，建立了能真实地反映岩石的强度参数内摩擦角和粘聚力与内变量的关系，即，修正的cwfs模型。在动态荷载作用下，岩石的强度和刚度往往会随应变率的增加而增加，因此，可用强度或刚度动态增强因子模型来描述岩石强度或刚度与应变率的关系。
[0089]
如图1所示，本发明一种适用于地震荷载下岩石动态本构模型的构建方法具体包括以下步骤：
[0090]
步骤一：基于次加载面理论和drucker
‑
prager屈服准则，构建岩石次加载面的函数以及岩石相似中心面的函数；
[0091]
步骤二：基于试验数据绘制岩石在不同围压下的偏应力
‑
轴向应变关系曲线、不同应变率下岩石单轴压缩的应力应变曲线以及岩石循环加卸载试验应力
‑
应变曲线；
[0092]
步骤三：将修正cwfs模型与动态增强因子模型进行耦合，以获取岩石的粘聚力的耦合计算模型，根据该耦合计算模型以及基于岩石弹性模量的损伤和应变率效应的耦合函数，并结合偏应力
‑
轴向应变关系曲线、应力应变曲线以及岩石循环加卸载试验应力
‑
应变曲线构建和求解弹塑性矩阵。
[0093]
步骤四：根据所述弹塑性矩阵构建岩石动态本构模型，同时更新步骤一中岩石正常屈服面与次加载面相似比的函数式，重复步骤一至步骤四，直至所构建的岩石动态本构模型计算的应力满足次加载面上的屈服函数。
[0094]
步骤五：基于有限元软件，对不同加载速率循环荷载试验进行数值模拟，验证此动态本构模型的适用性和准确性。
[0095]
具体而言，本发明中，可按如下步骤求解弹塑性矩阵：
[0096]
(1.1)求解屈服函数f和第一材料参数q(κ)的表达式：选择drucker
‑
prager屈服面作为次加载面理论的正常屈服面，则正常屈服面上的屈服函数可以表示为：
[0097][0098]
其中，β(κ)为第二材料参数，κ为内变量，i1和j2则分别是应力σ的第一不变量和应力的偏应力第二不变量。第二材料参数β(κ)和第一材料参数的表达式为：
[0099][0100][0101]
[0102][0103][0104][0105]
式中，和coh分别为岩石的内摩擦角和粘聚力，是随内变量κ变化的函数，采用了修正的cwfs模型，coh是随内变量κ和应变率变化的函数，采用了修正的cwfs模型和动态增强因子模型的耦合模型。c0和c
r
分别是粘聚力初始值和最终值，和分别是内摩擦角的初始值和最终值，和分别是粘聚力和内摩擦角开始变化时内变量的值，和分别是粘聚力和内摩擦角达到最终值时内变量的值，dε
p
是塑性应变增量，p为围压，f
c
为单轴抗压强度，a1和a2分别为第三材料参数和第四材料参数，tr为矩阵的迹，为静态荷载下的应变率，通常取1x10
‑5/s，d为强度动态增强因子模型的第五材料参数。
[0106]
单轴抗压强度f
c
可以根据岩石单轴压缩试验求解。根据三轴分级循环加卸载试验得出不同围压p下的应力应变曲线，把峰值前和峰值后的曲线分开考虑，在峰值前阶段，首先根据公式(6)计算不同围压p下从初始屈服到峰值间的塑性应变ε
p
，假设内变量κ在初始屈服时为0，在达到峰值时为1，然后得出不同围压p下内变量增量dκ与围压p的关系，进而求解出n的关系式；然后将内变量κ等间隔取值，将试验得到的偏应力根据对应的内变量κ的数值进行线性插值，从而得到不同围压p下相同内变量κ对应的偏应力，根据摩尔圆即可求出不同内变量κ对应的强度参数粘聚力coh和内摩擦角同样地，也可以求出峰后阶段求出不同内变量κ对应的强度参数粘聚力coh和内摩擦角根据峰值前的塑性应变ε
p
与峰后的塑性应变ε
p
的比值，即可把峰值前和峰后的内变量κ进行统一，从而得出强度参数粘聚力coh和内摩擦角与内变量κ的关系，进而可以求解c0，c
r
，，a1，a2。根据不同应变率下岩石单轴压缩的应力应变曲线，获取不同应变率下的岩石的峰值强度，通过强度增强因子模型拟合可以得出参数d。
[0107]
由于次加载面与正常屈服面保持几何相似的，因此次加载面上的函数为：
[0108][0109]
其中，为次加载面上应力，和则分别是次加载面上应力的第一不变量和次加载面上应力的偏应力第二不变量，r为次加载面与正常屈服面的大小之比，称为相似比，次加载面上应力由应力σ获得，由应力σ获得，为次加载面上的背应力。同理，相似中心面的函数为：
[0110][0111]
其中，s为相似中心，和则分别是相似中心面上相似中心s的第一不变量和相似中心面上相似中心s的偏应力第二不变量，r
s
为相似中心面与正常屈服面之比，也称相似中心比。
[0112]
(1.2)求解正常屈服面上背应力α的增量dα：对于岩土材料，循环荷载与动力问题应采用随动硬化或混合硬化，正常屈服面上背应力α采用非线性随动硬化准则，即：
[0113][0114]
其中，i为第六材料参数，r为第七材料参数，||||表示向量的模，dε
p
为应力σ对应的塑性应变增量。第六材料参数i和第七材料参数r可以通过岩石循环加卸载试验得出的应力
‑
轴向应变关系曲线在初始和反向加载阶段可以求出。
[0115]
(1.3)根据以下公式获得应力σ与相似中心s之矢量差次加载面上考虑正常屈服面上背应力α的应力以及相似中心s与正常屈服面上背应力α之矢量差
[0116][0117][0118][0119]
式中，σ
y
是应力σ在正常屈服面上的对偶应力，为相似中心s与次加载面上背应力之矢量差。
[0120]
(1.4)求解(1.4)求解和其表达式为：
[0121][0122][0123][0124]
[0125][0126][0127][0128][0129][0130][0131][0132]
式中，为次加载面上应力的平均值，s
m
为相似中心s的平均值。
[0133]
(1.5)求解弹塑性矩阵d
ep
：综合步骤(1.1)～(1.5)，可求得弹塑性矩阵d
ep
。
[0134][0135]
[0136][0137][0138][0139][0140][0141]
式中，“·”为点乘，d
el
为弹性矩阵，ν为泊松比，χ为相似中心比r
s
的最大值，χ值不超过1。δ是二阶单位矩阵，i是对称的四阶单位矩阵。m和u分别为第八材料参数和第九材料参数。a0和b分别为第十材料参数和第十一材料参数。是弹性模量，是随内变量κ和应变率变化的函数，更进一步而言，本发明中，将岩石弹性模量的损伤和应变率效应进行了耦合，从而得到的耦合函数。e
s
和e
e
分别是粘聚力初始值和最终值。κ
e
是弹性模量达到最终值时内变量的值。
[0142]
基于试验数据绘制岩石在不同围压下的偏应力
‑
轴向应变关系曲线，根据初始阶段偏应力
‑
轴向应变关系曲线直线段，可计算得出弹模初始值e
s
、泊松比ν，根据残余阶段偏应力
‑
轴向应变关系曲线直线段，可计算得出弹模最终值e
e
。同理，依据步骤(1.1)中的方法可以求解弹性模量达到最终值时内变量的值κ
e
。根据不同应变率下岩石单轴压缩的应力应变曲线，获取不同应变率下的岩石的弹模初始值，通过刚度增强因子模型拟合可以得出参数a0和b。u是控制应力点向正常屈服状态靠近的速率的参数，由屈服状态下岩石循环加卸载试验应力
‑
应变曲线的斜率初步确定，即弹塑性过渡状态。m的大小则影响滞回圈的宽度。调整u和m，直到地震荷载下岩石材料的动态本构模型较好地拟合岩石循环加卸载试验应力
‑
应变曲线为止。
[0143]
(1.6)求解相似比r1的具体函数式：
[0144]
相似比r1为：
[0145][0146]
其中：
[0147][0148][0149][0150]
式中，和分别为和的平均值，和分别为和的偏应力。将相似比r更新为相似比r1并返回步骤(1.1)，直至所构建的岩石动态本构模型计算的应力满足次加载面上的屈服函数，得到满足要求的动态本构模型。
[0151]
在所述步骤四中，可按如下步骤建立适用于地震荷载下岩石的动态本构模型的方程为：
[0152][0153]
式中，dσ为应力增量，dε为应变增量。
[0154]
在所述步骤五中，通过有限元商业软件(如abaqus)，把地震荷载下岩石材料的动态本构模型进行植入。对t
2y
6大理岩试样开展常规三轴压缩试验，将t
2y
6大理岩原状岩块切割打磨成上下端面水平且高度与直径比为2:1的标准圆柱样，本发明中岩样侧面光滑，直径为50mm，高度为100mm。岩样安装于试验仪器压力室内，利用伺服系统按静水压力条件对岩样施加围压至预定值并保持稳定，采用轴向位移控制加载，以一定速率对岩样施加轴向偏应力至岩样破坏，整个加载过程中，伺服测试系统自动记录施加的轴向偏应力值、轴向应变和径向应变。基于试验数据绘制t
2y
6大理岩在不同围压下的轴向偏应力
‑
轴向应变关系曲线、偏应力
‑
径向应变关系曲线，如图2所示，获得v和f
c
(单位为mpa)。
[0155]
进一步，对t
2y
6大理岩试样开展三轴分级循环加卸载试验，对标准圆柱样采用轴向位移控制加载，循环加、卸载的加载阶段的位移速率为0.06mm/min，其对应的轴向应力加载速率随围压的不同而稍有差别，范围为25～35mpa/min，加卸载试验的卸载段的轴向应力卸载速率为26mpa/min，卸载至偏压5mpa左右。整个加载过程中，伺服测试系统自动记录施加的轴向偏应力值、轴向应变和径向应变。处理后得出岩样的轴向偏应力
‑
轴向应变关系曲线和轴向偏应力
‑
径向应变关系曲线(包络线)，如图3、图4以及图5所示，经过处理，可以得出强度参数粘聚力coh和内摩擦角与内变量κ的关系如图7所示，同样也可以得出弹性模量e与内变量κ的关系如图9所示。可得出地震荷载下t
2y
6大理岩动态本构模型参数如表1所示，获得c0(单位为mpa)，c
r
(单位为mpa)，(单位为
°
)，(单位为
°
)，e
s
(单位为gpa)，e
e
(单位为gpa)，κ
e
，i，r、u、m、χ。
[0156]
此外，对t
2y
6大理岩试样开展单轴不同应变率动态压缩试验，对标准圆柱样采用轴向力控制加载，整个加载过程中，力随时间的增量都保持不变，即加载速率(或应变率)恒定，伺服测试系统自动记录施加的轴向应力值、轴向应变和径向应变。可以得出不同应变率
下弹性模量以及峰值强度，经过处理，可以得出峰值强度和弹性模量e与应变率的关系如图8和10所示，可得出地震荷载下t
2y
6大理岩动态本构模型参数如表1所示，获得a0、b、d。
[0157]
将e
s
(单位为gpa)，e
e
(单位为gpa)，v，f
c
(单位为mpa)，c0(单位为mpa)，c
r
(单位为mpa)，(单位为
°
)，(单位为
°
)，κ
e
，i，r、u、m、χ代入地震荷载下岩石材料的动态本构模型的方程，并植入限元商业软件(如abaqus)，对t
2y
6大理岩岩进行不同应变率循环加卸载数值模拟，包括静态应变率和应变率为4.3x10
‑5/s，通过监测点记录施加的轴向偏应力值、轴向应变和径向应变，与试验得出岩样在不同应变率下的轴向偏应力
‑
轴向应变关系曲线和轴向偏应力
‑
径向应变关系曲线进行对比验证，如图6和图11所示，数值模拟结果与试验结果呈现规律具有一致性，说明地震荷载下岩石材料的动态本构模型能准确地对t
2y
6大理岩不同应变率循环加载试验观察到的基本力学特性予以描述。由此认为本发明提出的岩石在地震荷载下的动态本构模型的方法结果准确，且力学意义明确。
[0158]
表1 t
2y
6大理岩参数表
[0159][0160]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

技术特征：
1.一种适用于地震荷载下岩石动态本构模型的构建方法，其特征在于，包括以下步骤：s1基于次加载面理论和drucker
‑
prager屈服准则，构建岩石次加载面的函数以及岩石相似中心面的函数；s2基于试验数据绘制岩石在不同围压下的偏应力
‑
轴向应变关系曲线、不同应变率下岩石单轴压缩的应力应变曲线以及岩石循环加卸载试验应力
‑
应变曲线；s3将修正cwfs模型与动态增强因子模型进行耦合，以获取岩石的粘聚力的耦合计算模型，根据该耦合计算模型以及基于岩石弹性模量的损伤和应变率效应的耦合函数，并结合偏应力
‑
轴向应变关系曲线、应力应变曲线以及岩石循环加卸载试验应力
‑
应变曲线构建和求解弹塑性矩阵；s4根据所述弹塑性矩阵构建岩石动态本构模型，同时更新步骤s1中岩石正常屈服面与次加载面相似比的函数式，重复步骤s1至步骤s4，直至所构建的岩石动态本构模型计算的应力满足次加载面上的屈服函数。2.根据权利要求1所述的方法，其特征在于，步骤s1具体包括以下步骤：s11选择drucker
‑
prager屈服面作为次加载面理论的正常屈服面，构建正常屈服面上的屈服函数的屈服函数其中，β(κ)为第二材料参数，κ为内变量，i1和j2则分别是应力σ的第一不变量和应力的偏应力第二不变量；s12基于次加载面理论，构建岩石次加载面的函数s12基于次加载面理论，构建岩石次加载面的函数其中，为次加载面上的函数，为次加载面上应力，和则分别是次加载面上应力的第一不变量和次加载面上应力的偏应力第二不变量，r为次加载面与正常屈服面的相似比，β(κ)为第二材料参数，κ为内变量，为第一材料参数；s13构建岩石相似中心面的函数s13构建岩石相似中心面的函数其中，s为相似中心，和则分别是相似中心面上相似中心s的第一不变量和相似中心面上相似中心s的偏应力第二不变量，r
s
为相似中心面与正常屈服面之比。3.根据权利要求2所述的方法，其特征在于，步骤s11中，所述第二材料参数β(κ)的表达式为：所述第一材料参数的表达式为：
式中，coh为岩石的粘聚力，该粘聚力是随内变量κ和应变率变化的函数。4.根据权利要求1所述的方法，其特征在于，步骤s3中，所述岩石的粘聚力的耦合计算模型为：式中，为岩石的内摩擦角，该内摩擦角是随内变量κ变化的函数，其具体计算模型为：所述内变量κ的具体计算模型为：所述内变量κ的具体计算模型为：其中，c0为粘聚力初始值，c
r
为粘聚力最终值，为内摩擦角的初始值，为内摩擦角的最终值，是粘聚力开始变化时内变量的值，是内摩擦角开始变化时内变量的值，是粘聚力达到最终值时内变量的值，是内摩擦角达到最终值时内变量的值，dε
p
是塑性应变增量，p是围压，f
c
是单轴抗压强度，a1和a2分别为第三材料参数和第四材料参数，tr为矩阵的迹，为静态荷载下的应变率，d为强度动态增强因子模型的第五材料参数，n是围压p的函数，其用于表示内变量k的增长速率。5.根据权利要求4所述的方法，其特征在于，步骤s3中，所述单轴抗压强度f
c
根据岩石单轴压缩试验求解，根据三轴分级循环加卸载试验得出不同围压p下的应力应变曲线，把峰值前和峰值后的曲线分开考虑，在峰值前阶段，首先计算不同围压p下从初始屈服到峰值间的塑性应变ε
p
，假设内变量κ在初始屈服时为0，在达到峰值时为1，然后得出不同围压p下内变量增量dκ与围压p的关系，进而求解出n的关系式；然后将内变量κ等间隔取值，将试验得到的偏应力根据对应的内变量κ的数值进行线性插值，从而得到不同围压p下相同内变量κ对应的偏应力，根据摩尔圆即可求出不同内变量κ对应的强度参数粘聚力coh和内摩擦角
接着，求出峰后阶段求出不同内变量κ对应的强度参数粘聚力coh和内摩擦角根据峰值前的塑性应变ε
p
与峰后的塑性应变ε
p
的比值，把峰值前和峰后的内变量κ进行统一，从而得出强度参数粘聚力coh和内摩擦角与内变量κ的关系，进而可以求解c0，c
r
，a1，a2，根据不同应变率下岩石单轴压缩的应力应变曲线，获取不同应变率下的岩石的峰值强度，通过强度增强因子模型拟合得出参数d。6.根据权利要求4所述的方法，其特征在于，所述步骤s3具体包括以下步骤：s31求解正常屈服面上背应力α的增量dα：式中，i为第六材料参数，r为第七材料参数，|| ||表示向量的模，dε
p
为应力σ对应的塑性应变增量；s32获取应力σ与相似中心s的矢量差次加载面上考虑正常屈服面上背应力α的应力以及相似中心s与正常屈服面上背应力α的矢量差以及相似中心s与正常屈服面上背应力α的矢量差以及相似中心s与正常屈服面上背应力α的矢量差以及相似中心s与正常屈服面上背应力α的矢量差式中，σ
y
是应力σ在正常屈服面上的对偶应力，为相似中心s与次加载面上背应力的矢量差；s33构建并求解弹塑性矩阵d
ep
，该弹塑性矩阵d
ep
的计算模型如下：式中，式中，式中，
其中，“·”为点乘，d
el
为弹性矩阵，ν为泊松比，χ为相似中心比r
s
的最大值，δ是二阶单位矩阵，i是对称的四阶单位矩阵，m和u分别为第八材料参数和第九材料参数，a0和b分别为第十材料参数和第十一材料参数，是弹性模量，是随内变量κ和应变率变化的函数，e
s
和e
e
分别是粘聚力初始值和最终值，κ
e
是弹性模量达到最终值时内变量的值。7.根据权利要求5所述的方法，其特征在于，步骤s3中，基于试验数据绘制岩石在不同围压下的偏应力
‑
轴向应变关系曲线，根据初始阶段偏应力
‑
轴向应变关系曲线直线段，计算得出弹模初始值e
s
、泊松比ν，根据残余阶段偏应力
‑
轴向应变关系曲线直线段，计算得出弹模最终值e
e
；根据不同应变率下岩石单轴压缩的应力应变曲线，获取不同应变率下的岩石的弹模初始值，通过刚度增强因子模型拟合得出中的参数a0和b；u为控制应力点向正常屈服状态靠近的速率的参数，通过屈服状态下岩石循环加卸载试验应力
‑
应变曲线的斜率确定。8.根据权利要求1所述的方法，其特征在于，所述岩石正常屈服面与次加载面相似比的函数式中，将更新后的次加载面与正常屈服面的相似比r更新为r1，所述相似比r1的计算公式如下：其中，其中，
式中，为的平均值，为的平均值，为的偏应力，为的偏应力，β(κ)为第二材料参数，κ为内变量，为第一材料参数，为应力σ与相似中心s的矢量差，为相似中心s与正常屈服面上背应力α的矢量差。9.根据权利要求1所述的方法，其特征在于，所述岩石动态本构模型的计算公式如下：其中，dσ为应力增量，dε为应变增量，为次加载面上的函数。10.根据权利要求1所述的方法，其特征在于，还包括以下步骤：基于有限元软件，对不同加载速率循环荷载试验进行数值模拟，验证所述岩石动态本构模型的适用性和准确性。
技术总结
本发明属于岩体工程技术领域，并具体公开了一种适用于地震荷载下岩石动态本构模型的构建方法。包括：基于次加载面理论和Drucker


技术研发人员：周永强 盛谦 付晓东
受保护的技术使用者：中国科学院武汉岩土力学研究所
技术研发日：2021.02.05
技术公布日：2021/6/25