CN110174307B - 基于横观各向同性岩体三维蠕变模型的参数辨识方法 - Google Patents

基于横观各向同性岩体三维蠕变模型的参数辨识方法 Download PDF

Info

Publication number
CN110174307B
CN110174307B CN201910338061.1A CN201910338061A CN110174307B CN 110174307 B CN110174307 B CN 110174307B CN 201910338061 A CN201910338061 A CN 201910338061A CN 110174307 B CN110174307 B CN 110174307B
Authority
CN
China
Prior art keywords
creep
bedding
formula
rock sample
stress
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
Application number
CN201910338061.1A
Other languages
English (en)
Other versions
CN110174307A (zh
Inventor
李列列
管俊峰
姚贤华
何双华
高勇伟
高新亮
张玉龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
North China University of Water Resources and Electric Power
Original Assignee
North China University of Water Resources and Electric Power
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by North China University of Water Resources and Electric Power filed Critical North China University of Water Resources and Electric Power
Priority to CN201910338061.1A priority Critical patent/CN110174307B/zh
Publication of CN110174307A publication Critical patent/CN110174307A/zh
Application granted granted Critical
Publication of CN110174307B publication Critical patent/CN110174307B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N3/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N3/08—Investigating strength properties of solid materials by application of mechanical stress by applying steady tensile or compressive forces
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/0014—Type of force applied
    • G01N2203/0016—Tensile or compressive
    • G01N2203/0019—Compressive
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/0058—Kind of property studied
    • G01N2203/0069—Fatigue, creep, strain-stress relations or elastic constants
    • G01N2203/0071—Creep
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/02—Details not specific for a particular testing method
    • G01N2203/0202—Control of the test
    • G01N2203/0212—Theories, calculations
    • G01N2203/0218—Calculations based on experimental data
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/02—Details not specific for a particular testing method
    • G01N2203/025—Geometry of the test
    • G01N2203/0252—Monoaxial, i.e. the forces being applied along a single axis of the specimen
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/02—Details not specific for a particular testing method
    • G01N2203/025—Geometry of the test
    • G01N2203/0256—Triaxial, i.e. the forces being applied along three normal axes of the specimen
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06—Power analysis or power optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Analytical Chemistry (AREA)
  • Computer Hardware Design (AREA)
  • Immunology (AREA)
  • Computing Systems (AREA)
  • General Health & Medical Sciences (AREA)
  • Biochemistry (AREA)
  • Evolutionary Computation (AREA)
  • Pathology (AREA)
  • Chemical & Material Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Geometry (AREA)
  • Operations Research (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明公开了一种基于横观各向同性岩体三维蠕变模型的参数辨识方法，旨在解决现有技术中存在的岩石材料的各向异性蠕变本构模型难以全面反映出横观各向同性岩体的蠕变力学特性的技术问题。本发明首次提出了横观各向同性岩体三维模型蠕变参数的辨识方法，可实现采用一套蠕变参数同时表现出横观各向同性岩体不同层理方向蠕变特征的目的，应用范围广泛。

Description

基于横观各向同性岩体三维蠕变模型的参数辨识方法
技术领域
本发明涉及岩石力学与工程技术领域，具体涉及一种基于横观各向同性岩体三维蠕变模型的参数辨识方法。
背景技术
层状岩体在地壳中广泛存在，其在力学上具有横观各向同性的特点，岩体的破坏和失稳与其内部的损伤演化有着密切的关系。
随着地下工程的开发建设，层状岩体的力学性质（包括蠕变特性）对隧道、地下厂房、水电等工程围岩变形和应力分布，尤其是长期稳定性有重要的影响。由于层状岩体特殊的构造性质，其力学性质与各向同性岩石相比具有较大的差异，如随着层理倾角的不同，其蠕变特性具有显著的各向异性。如溪洛渡工程坝肩两岸分别存在N200W和近EW向两组基体结构面，这两组结构面在岩体中大量发育，构成岩体的基本层次，两岸岩体沿结构面发生滑移错动挤压坝体，对坝体安全性构成严重威胁；锦屏一级水电站坝址区地层岩性、层厚变化大，岩体组合复杂，其中大理岩夹杂工程地质特性较差的绿泥石组成的层状岩体，主要分布在两岸坝肩、水垫塘及地下厂房等部位，影响到这些工程部位的稳定性，尤其对坝区右岸而言，大理岩中呈透镜状、团块状分布的绿片岩，其片理发育且倾向坡外，形成了层间弱面，成为影响右岸抗滑稳定的主要因素；内蒙古平庄西露天矿某处边坡，其下部的缓倾角砂岩—泥岩软弱夹层在降水的长期浸泡下，岩体软化，沿软弱面发生滑移，造成了滑坡事故。由此可以看出，岩体蠕变特性对工程建筑物的安全有着重要影响。从理论上分析评价这些建筑物的安全性是非常必要的，而构建科学、合理的蠕变本构方程是进行理论分析的基础。这其中就需要选用合理的蠕变参数进行分析预测，并利用分析结果指导工程设计施工，达到确保工程结构的安全性的目的。
现有的针对岩石材料的各向异性蠕变本构模型多是基于塑性理论，其更适用于金属等发生塑性破坏的材料。众所周知，岩体的破坏展示出准脆性性质，而非完全塑性。而现有模型研究成果多根据试验现象，通过构建一维蠕变本构模型，将不同方向蠕变参数简单辨识为不同的相互独立的值，而并未建立真正意义上的横观各向同性三维蠕变本构模型。
此外，由于岩体层理弱面的存在，现有的各向同性蠕变模型难以全面反映出横观各向同性岩体的蠕变力学特性。
因此，亟需研究建立适用于可描述岩体准脆性断裂破坏特性的各向异性黏弹塑性元件模型，进而辨识和解算出能够更为真实地反映岩体断裂破坏本质的相关表征参数，以应用于各类岩体工程建设的安全性预测评价，为其设计及决策提供更加精准可靠的科学依据。
发明内容
本发明要解决的技术问题是提供一种基于横观各向同性岩体三维蠕变模型的参数辨识方法，以期解决现有技术中所存在的岩石材料的各向异性蠕变本构模型难以全面反映出横观各向同性岩体的蠕变力学特性的技术难题。
为解决上述技术问题，本发明采用了如下的技术思路（参见图1）：
横观各向同性岩体由于层理弱面的存在，表现出各向异性的特点。因此，不同于各向同性岩石，在确定横观各向同性岩体五个弹性常数时，根据弹性本构方程的特点，至少需要根据平行组、垂直组两个试件的单轴或三轴压缩实验结果来计算弹性常数。蠕变参数也是如此，依据横观各向同性岩体三维蠕变本构方程的特点，根据平行组、垂直组两个试件的蠕变试验结果，首次分析得出横观各向同性岩体的蠕变参数。
具体采用的技术方案如下：
设计一种基于横观各向同性岩体三维蠕变模型的参数辨识方法，包括以下步骤：
（1）将岩石试样制成层理水平岩样和层理垂直岩样，进行横观各向同性岩体的三轴或单轴压缩蠕变试验；
（2）采用陈氏加载法进行处理步骤（1）的试验结果，得层理水平岩样和层理垂直岩样的分级蠕变试验曲线；
（3）三维空间中，设xoy为层理面所在平面，采用以下层状岩体的三维蠕变本构方程，进行相应应变计算：
Figure 359778DEST_PATH_IMAGE001
——式（I）；
Figure 166191DEST_PATH_IMAGE002
——式（II）；
Figure 920521DEST_PATH_IMAGE003
——式（III）；
Figure 657532DEST_PATH_IMAGE004
——式（IV）；
Figure 243235DEST_PATH_IMAGE005
——式（V）；
Figure 486128DEST_PATH_IMAGE006
——式（VI）；
式（I）至式（VI）中：ε x 为x方向的正应变；ε y 为y方向的正应变；ε z 为z方向的正应变；γ yz 为yoz平面的剪应变；γ xz 为xoz平面的剪应变；γ xy 为xoy平面的剪应变；μ hh 、μ hv 分别为平行和垂直层理方向的泊松比；Εh M、ηh M、Εh K、ηh K，及对应的n 1 、n 2 、n 3 、n 4 分别为垂直层理面方向的八个相互独立的拉压蠕变参数；t为时间；σ x 为方向的正应力，σ y 为方向的正应力，σ z 为方向的正应力；
（4）对于层理水平岩样，xoy平面为各向同性平面，则σ x= σ y= σ 3 ，σ z= σ 1 ，ε z= ε0轴，其应变采用式（III）计算，式（III）变换写为式（VII）：
Figure 727754DEST_PATH_IMAGE007
——式（VII）；
对于层理垂直岩样，σ x= σ y= σ 3 ，σ x= σ 1 ，ε x= ε90轴，其应变采用式（I）计算，式（I）变换写为式（VIII）：
Figure 2877DEST_PATH_IMAGE008
——式（VIII）；
式中，σ 1 为轴向大主应力；σ 3 为环向小主应力；
（5）基于式（VII），结合层理水平岩样的蠕变试验曲线，通过Levenberg-Marquardt非线性最小二乘法进行拟合得到参数组合Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 ；
（6）若层理水平岩样和垂直岩样采用相同的应力水平进行蠕变试验，即层理水平岩样应力水平σ H 等于垂直岩样的应力水平σ V 则将所述Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 代入式（VIII）中；根据层理垂直岩样的蠕变试验曲线结果，拟合得到Εh M、ηh M、Εh K和ηh K；
（7）将所述Εh Mn 1 、ηh Mn 2 、Εh Kn 3 、ηh Kn 4 和所述Εh M、ηh M、Εh K、ηh K的计算结果相除，得出n 1 、n 2 、n 3 和n 4 ，得横观各向同性岩体的八个独立的拉压蠕变参数。
优选的，在所述步骤（6）中，若层理水平岩样和垂直岩样采用不同的应力水平进行蠕变试验，即σ H ≠σ V ，需根据层理水平岩样的试验结果，通过回归分析似计算得到层理垂直岩样应力水平对应的Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 。
进一步的，为避免出现过大值、过小值甚至负值使结果失真，所述回归分析中回归曲线不选用线性和/或多项式方程。
优选的，所述回归分析中回归曲线采用非线性方程。
优选的，所述非线性方程为指数方程、幂函数方程、对数方程中的任一种。
与现有技术相比，本发明的有益技术效果在于：
1. 本发明首次提出了基于横观各向同性岩体三维模型蠕变参数的辨识方法。
2. 本发明所依据的模型综合考虑了平行和垂直层理面方向蠕变力学行为的差异性，依据横观各向同性岩体的三维蠕变本构模型，能够辨识并解算出八个独立的拉压蠕变参数。
3. 本发明蠕变参数识别方法可实现采用一套蠕变参数同时表征出横观各向同性岩体不同层理方向的蠕变特征的目的。
4. 采用本发明所提方法，既可采用传统的蠕变试验设计方案（即根据不同层理倾角千枚岩荷载峰值，分级加载）进行不同应力水平的蠕变试验；也可采用采用相同的应力水平进行蠕变试验，应用范围更为广泛。
5. 本发明建立的蠕变参数识别方法通过选用合理的蠕变参数进行分析，并利用分析结果指导工程设计施工，可确保工程结构的安全性（如对四川杜家山千枚岩隧道进行稳定性分析），还可对工程的长期安全性进行预测评价，对岩石力学理论以及地下工程建设具有重大的应用价值。
附图说明
图1为横观各向同性岩体三维模型蠕变参数识别流程图；
图2为蠕变模型示意图；
图3为坐标系示意图；
图4为千枚岩层理水平岩样的轴向蠕变分级曲线图；
图5为千枚岩层理垂直岩样的轴向蠕变分级曲线图；
图6为千枚岩蠕变参数回归分析图；
图7为千枚岩层理水平岩样的蠕变试验及发展模型的拟合曲线图；
图8为千枚岩层理垂直岩样的蠕变试验及发展模型的拟合曲线图；
图9为云母片岩层理水平岩样的轴向蠕变分级曲线图；
图10为云母片岩层理垂直岩样的轴向蠕变分级曲线图；
图11为云母片岩蠕变参数回归分析图；
图12为云母片岩层理水平岩样的蠕变试验及发展模型的拟合曲线图；
图13为云母片岩层理垂直岩样的蠕变试验及发展模型的拟合曲线图。
具体实施方式
下面结合附图和实施例来说明本发明的具体实施方式，但以下实施例只是用来详细说明本发明，并不以任何方式限制本发明的范围。
在以下实施例中所涉及的材料如无特别说明，均为市售常规材料；所涉及的步骤或试验方法，如无特别说明，均为常规方法；所涉及的名称参数，如无特别说明，均为本领域常规术语。
实施例一：层状岩体的三维蠕变本构方程的构建
1. 层状岩体的蠕变特性具有以下特征：
（1）加载初期具有瞬时弹性变形，模型中应含有弹性元件；
（2）随着时间的增加，蠕变速率逐渐降低，并逐渐近似趋近于某个常数，故模型中应具有粘性元件；
（3）在低应力水平下，蠕变曲线达到某个极限值后，其应变值保持不变，因此模型中应包含弹性元件与粘性元件的并联；
可得描述层岩岩体的蠕变元件模型如图2所示。
2. 对图2元件模型进行力学分析。
（1）Maxwell体中弹性元件与粘性元件串联，具有以下关系：
Figure DEST_PATH_IMAGE009
——式（a）；
Figure 708665DEST_PATH_IMAGE010
——式（b）；
式中，D为算子，D=d/dt。
（2）Kelvin体中弹性元件与粘性元件并联，具有以下关系：
Figure DEST_PATH_IMAGE011
——式（c）；
Figure 853951DEST_PATH_IMAGE012
——式（d）；
式（a）～式（d）中，E M ，E K 分别为Maxwell体和Kelvin体的弹性参数；η M ，η K 分别为Maxwell体和Kelvin体的粘性参数。
图2中的一维蠕变模型由Maxwell体和Kelvin体串联而成，可得：
Figure DEST_PATH_IMAGE013
——式（e）；
将式（b）、（d）代入式（e），得：
Figure 582873DEST_PATH_IMAGE014
——式（f）；
计蠕变柔量：
Figure DEST_PATH_IMAGE015
，则式（f）可写为：
Figure 724004DEST_PATH_IMAGE016
——式（g）；
通过对式（f）中的算子D进行积分，可得到最终的一维蠕变本构方程：
Figure DEST_PATH_IMAGE017
——式（h）；
3. 各向同性材料蠕变本构方程
根据一维条件下建立的蠕变本构方程，通过常泊松比方法将其推广到三维应力状态。
思路为：物体蠕变过程中假设泊松比不随时间和应力改变，并等于弹性阶段的泊松比，μ(σ,t)=μ为不变的实常数。基于此假定将蠕变方程从一维状态推广到三维状态。
各向同性材料在三维应力空间的应力应变关系为：
Figure 35031DEST_PATH_IMAGE018
——式（i）；
式中：{ε e }为弹性应变矩阵，{σ e }为三维应力矩阵，{H}为柔度矩阵，分别可表示为：
Figure DEST_PATH_IMAGE019
Figure 806678DEST_PATH_IMAGE020
——式（j）；
矩阵（j）中E、μ分别为各向同性材料的弹性模量和泊松比。
式（i）也可写为：
Figure DEST_PATH_IMAGE021
——式（k）；
式中：
Figure 819633DEST_PATH_IMAGE022
。
保持泊松比矩阵[U]不发生变化，采用微分算子替换的方法，将蠕变模型中的柔度模量J(t)替换线弹性本构，即式（k）中的柔度系数1/E，得到蠕变本构方程的基本式：
Figure DEST_PATH_IMAGE023
——式（l）；
根据前文的推导
Figure 515188DEST_PATH_IMAGE024
，代入式（l），得
Figure DEST_PATH_IMAGE025
——式（m）；
将式（m）采用张量形式表达，可得：
Figure 867672DEST_PATH_IMAGE026
——式（n）；
式中：i=x、y、z，j=x、y、z，σm=1/3(σx+σy+σz)，δij为Kronec ker 符号。
4. 横观各向同性岩体蠕变本构方程
在横观各向同性条件下，假设横观各向同性岩体的自然坐标系为Oxyz，横观各向同性面所在平面为xOy平面，z轴为xOy平面法向方向，x、y、z服从右手法则，如图3所示。
基于图3，可得三维应力状态下横观各向同性岩体的弹性应力应变本构关系为：
Figure DEST_PATH_IMAGE027
——式（o）；
式中：
Figure 872537DEST_PATH_IMAGE028
其中，E h 、E v 、μ hh 、μ hv 、G’为5个弹性常数参数，E h 、E v 分别为平行和垂直层理方向的弹性模量，μ hh 、μ hv 分别为平行和垂直层理方向的泊松比，G’为垂直层理方向的剪切模量，{ε}={ε x 、ε y 、ε z 、γ yz 、γ xz 、γ xy } T ，{σ}={σ x 、σ y 、σ z 、τ yz 、τ xz 、τ xy } T 。
横观各向同性岩体具有五个独立的弹性参数。Saint—Venant提出了近似弹性参数表达式（p），其为四个弹性参数。
Figure DEST_PATH_IMAGE029
——式（p）；
假设E v =nE h ，并结合式（p），则式（o）中的柔度矩阵[S]可写为：
Figure 123521DEST_PATH_IMAGE030
——式（q）；
式中：
Figure DEST_PATH_IMAGE031
针对横观各向同性岩体的特点，基于常泊松比思想，假定μhh和μhv不随应力和时间改变，均等于弹性阶段的常数。采用与各向同性材料推广到三维相同的思想，保持泊松比矩阵不发生变化，将式（l）中的各向同性材料的泊松比矩阵[U]采用式（q）所示的横观各向同性材料泊松比矩阵[Utra]代替，则可得层状岩体的三维蠕变本构方程如下：
Figure 544138DEST_PATH_IMAGE032
——式（r）。
实施例二：横观各向同性岩体三维蠕变模型的应用
基于上述模型，可实现采用一套蠕变参数同时表现出横观各向同性岩体不同层理方向的蠕变特征的目的，具体包括以下步骤：
（1）基于横观各向同性岩体的层理面与水平面的不同夹角，岩石试样分别加工制作成层理水平岩样和层理垂直岩样两类，进行横观各向同性岩体的三轴压缩蠕变试验；
（2）将层理水平岩样与层理垂直岩样的试验结果，采用陈氏加载法（参考文献：陈宗基，康文法．岩石的封闭应力、蠕变和扩容及本构方程）进行处理，进一步得到相应的层理水平岩样与层理垂直岩样的分级蠕变试验曲线；
（3）三维空间中，设xoy为层理面所在平面，将实施例一中的式（r）中蠕变柔量进行积分计算，并考虑平行和垂直层理方面蠕变特性的差异，横观各向同性材料三维蠕变模型由垂直层理方向的Εh M、ηh M、Εh K、ηh K，及对应的n 1 、n 2 、n 3 、n 4 共8个独立拉压蠕变参数，将式（r）展开，并采用类比的思想，可得横观各向同性材料的蠕变方程表达式，采用如下层状岩体的三维蠕变本构方程，进行相应应变计算：
Figure 79024DEST_PATH_IMAGE001
——式（I）；
Figure 192474DEST_PATH_IMAGE002
——式（II）；
Figure 117704DEST_PATH_IMAGE003
——式（III）；
Figure 152132DEST_PATH_IMAGE004
——式（IV）；
Figure 213629DEST_PATH_IMAGE005
——式（V）；
Figure 497980DEST_PATH_IMAGE006
——式（VI）；
式（I）至式（VI）中：ε x 为x方向的正应变；ε y 为y方向的正应变；ε z 为z方向的正应变；γ yz 为yoz平面的剪应变；γ xz 为xoz平面的剪应变；γ xy 为xoy平面的剪应变；μ hh 、μ hv 分别为平行和垂直层理方向的泊松比；Εh M、ηh M、Εh K、ηh K，及对应的n 1 、n 2 、n 3 、n 4 分别为垂直层理面方向的八个相互独立的拉压蠕变参数；t为时间。
（4）对于层理水平岩样，xoy平面为各向同性平面，则σ x= σ y= σ 3 ，σ z= σ 1 ，ε z= ε0轴，其应变采用式（III）计算，式（III）通过变换可写为式（VII）：
Figure 972823DEST_PATH_IMAGE007
——式（VII）；
式中，σ 1 为轴向大主应力；σ 3 为环向小主应力。
对于层理垂直岩样，σ x= σ y= σ 3 ，σ x= σ 1 ，ε x= ε90轴，其应变采用式（I）计算，式（I）通过变换可写为式（VIII）：
Figure 822DEST_PATH_IMAGE008
——式（VIII）；
（5）基于式（VII），结合层理水平岩样的蠕变试验曲线，通过Levenberg-Marquardt非线性最小二乘法进行拟合可分别得到参数组合Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 ；
（6）若层理水平岩样和垂直岩样采用相同的应力水平进行蠕变试验，即层理水平岩样应力水平σ H 等于垂直岩样的应力水平σ V ，则将第（5）步拟合得到的Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 结果代入式（VIII）中；再根据层理垂直岩样的蠕变试验曲线结果，拟合得到Εh M、ηh M、Εh K和ηh K；
（7）最后将Εh Mn 1 、ηh Mn 2 、Εh Kn 3 、ηh Kn 4 和对应Εh M、ηh M、Εh K、ηh K的计算结果，将二者相除，从而得出n 1 、n 2 、n 3 和n 4 。最终得到横观各向同性岩体的八个独立的拉压蠕变参数Εh M、ηh M、Εh K、ηh K、n 1 、n 2 、n 3 、n 4 。
试验例1：
取汶川县桃坪隧道进洞口的千枚岩作为试验材料，该岩石材料层状构造显著，为一种典型的横观各向同性岩体。采用实施例的方法验证所构建的横观各向同性岩体三维蠕变模型。
具体步骤如下：
（1）基于层理面与水平面的不同夹角，岩石试样被别加工制作成层理水平岩样和层理垂直岩样两类。
（2）进行三轴压缩蠕变试验：
采用5MPa的恒定围压，分级加载轴向压力的加载制度。其中，层理水平岩样轴向压力依次为9.4MPa、15.1MPa、20.8MPa、26.5MPa、32.2MPa；层理垂直岩样轴向压力依次为26.0MPa、35.6MPa、45.6MPa、55.5MPa、65.4MPa、75.1MPa、79.3MPa。
（3）将层理水平岩样与层理垂直岩样的试验结果采用陈氏加载法进行处理，室内试验为每间隔一段时间测量一个应变值，将不同的应变值连起来，进一步可得层理水平岩样与层理垂直岩样的分级蠕变试验曲线，如图4和图5所示。
由图4和图5可以看出，平行组和垂直组千枚岩在加载瞬间都会产生一定量的瞬时应变，此后在轴向应力作用下，产生随时间增长的蠕变变形。应力水平较低时，蠕变速率短时间内趋于零值；随着应力水平的提高，蠕变速率随时间的增长而衰减，减速蠕变时间增加，最后趋于一个恒定的非零值，如果蠕变时间足够长，蠕变变形达到一定值后最终使岩样发生破坏。
由此可见，千枚岩具有三轴应力作用下产生瞬时应变、衰减蠕变和稳定蠕变的典型蠕变行为特性。
（4）采用实施例中的蠕变参数识别方法与流程，基于步骤（3）所得层理水平岩样的蠕变试验结果，采用式（VII）进行拟合分析，可得出层理水平岩样在9.4MPa、15.1MPa、20.8MPa、26.5MPa、32.2MPa等不同应力水平下参数组合Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 的具体数据，确定结果如表1所示：
表1层理水平岩样蠕变参数
Figure DEST_PATH_IMAGE033
。
（5）根据层理水平岩样的试验结果，通过回归分析似计算得到层理垂直岩样应力水平对应的Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 。具体回归分析时，为避免出现过大值、过小值甚至负值而使结果失真，回归曲线一般情况下不选用线性和多项式方程，而采用指数方程、幂函数方程、对数方程等非线性方程。层理水平岩样的蠕变参数回归如图6所示。
由图6可以看出，Εh Mn 1 随着应力水平的提高有线性增加的趋势；ηh Mn 2 随着应力水平的提高呈指数函数下降；Εh Kn 3 和ηh Kn 4 随着应力水平呈非线性先升高后下降的趋势。
根据回归方程计算得到层理垂直岩样对应应力水平下的Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 的具体结果如表2所示。
表2 参数一览表
Figure 995454DEST_PATH_IMAGE034
。
（6）将表2中的分析结果代入式（VIII），并利用式（VIII）对层理垂直岩样的蠕变曲线进行拟合得到的蠕变参数如表3所示。
表3 参数一览表
Figure DEST_PATH_IMAGE035
。
（7）根据表2和表3中的参数预测结果，将二者相除，可计算得到n 1 、n 2 、n 3 和n 4 。 最终可得到横观各向同性岩体八个独立的蠕变参数，如表4所示。
表4 层理垂直岩样蠕变参数
Figure 185127DEST_PATH_IMAGE036
。
（8）将千枚岩蠕变试验所得值（注：图4、图5为试验值）与基于实施例的模型计算所得值拟合（注：拟合过程为采用与室内试验相同的时间t，即图中的横坐标值，代入蠕变方程式（VII）和式（VIII）中，方程中的蠕变参数如表4所示，从而对应每个时间t，计算出一个应变值，然后与试验值进行对比），结果如图7和图8所示。
由图7和图8可以看出，实施例二的三维蠕变损伤本构方程得出的计算值与试验值吻合良好，该模型可较好地反映了千枚岩的蠕变特性。证明了本发明横观各向同性岩体三维蠕变模型的表征有效、合理，能够反映出其真实的蠕变特性。
试验例2：
取国家重点工程丹巴水电站的云母片岩作为试验材料，该岩石材料层状构造显著，为一种典型的横观各向同性软岩体，对其蠕变损伤的分析直接影响丹巴水电站安全运营的理论分析结果。采用实施例的方法验证所构建的横观各向同性岩体三维蠕变模型。
具体步骤如下：
（1）基于层理面与水平面的不同夹角，岩石试样被别加工制作成层理水平岩样和层理垂直岩样两类。
（2）进行单轴压缩蠕变试验：
采用0MPa的围压（即无围压），分级加载轴向压力的加载制度。其中，层理水平岩样轴向压力依次为6.0MPa、10.0MPa、16.0MPa；层理垂直岩样轴向压力依次为3.0MPa、5.0MPa、7.0MPa、9.0MPa。
（3）将层理水平岩样与层理垂直岩样的试验结果采用陈氏加载法进行处理，室内试验为每间隔一段时间测量一个应变值，将不同的应变值连起来，进一步可得层理水平岩样与层理垂直岩样的分级蠕变试验曲线，如图9和图10所示。
由图9和图10可以看出，平行组和垂直组云母片岩在加载瞬间都会产生一定量的瞬时应变，此后在轴向应力作用下，产生随时间增长的蠕变变形。应力水平较低时，蠕变速率短时间内趋于零值；随着应力水平的提高，蠕变速率随时间的增长而衰减，减速蠕变时间增加，最后趋于一个恒定的非零值，如果蠕变时间足够长，蠕变变形达到一定值后最终使岩样发生破坏，应力越高，稳态蠕变速率越高。
由此可见，云母片岩在单轴应力作用下会产生瞬时应变、衰减蠕变和稳定蠕变的典型蠕变行为特性。
（4）采用实施例中的蠕变参数识别方法与流程，基于步骤（3）所得层理水平岩样的蠕变试验结果，采用式（VII）进行拟合分析，可得出层理水平岩样在6.0MPa、10.0MPa、16.0MPa等不同应力水平下参数组合Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 的具体数据，确定结果如表5所示：
表5层理水平岩样蠕变参数
Figure DEST_PATH_IMAGE037
。
（5）根据层理水平岩样的试验结果，通过回归分析似计算得到层理垂直岩样应力水平对应的Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 。具体回归分析时，为避免出现过大值、过小值甚至负值使结果失真，回归曲线一般情况下不选用线性和多项式方程，而采用指数方程、幂函数方程、对数方程等非线性方程。层理水平岩样的蠕变参数回归如图11所示。
由图11可以看出，Εh Mn 1 随着应力水平的提高有线性增加的趋势；ηh Mn 2 随着应力水平的提高呈指数函数先下降后提高的趋势；Εh Kn 3 和ηh Kn 4 随着应力水平呈指数函数先升高后下降的趋势。
根据回归方程计算得到层理垂直岩样对应应力水平下的Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 的具体结果如表6所示。
表6 参数一览表
Figure 147267DEST_PATH_IMAGE038
。
（6）将表6中的分析结果代入式（VIII），并利用式（VIII）对层理垂直岩样的蠕变曲线进行拟合得到的蠕变参数如表7所示。
表7 参数一览表
Figure DEST_PATH_IMAGE039
。
（7）根据表6和表7中的参数预测结果，将二者相除，可计算得到n 1 、n 2 、n 3 和n 4 。
最终可得到横观各向同性岩体八个独立的蠕变参数，如表8所示。
表8 层理垂直岩样蠕变参数
Figure 526427DEST_PATH_IMAGE040
。
（8）将云母片岩蠕变试验所得值（注：图9、图10为试验值）与基于实施例的模型计算所得值拟合（注：拟合过程为采用与室内试验相同的时间t，即图中的横坐标值，代入蠕变方程式（VII）和式（VIII）中，方程中的蠕变参数如表8所示，从而对应每个时间t，计算出一个应变值，然后与试验值进行对比），结果如图12和图13所示。
由图12和图13可以看出，实施例二的三维蠕变损伤本构方程得出的计算值与试验值吻合良好，该模型可较好地反映了云母片岩的蠕变特性。证明了本发明所构建的横观各向同性岩体三维蠕变模型的的表征有效、合理，能够反映出其真实的蠕变特性。
上面结合附图和实施例对本发明作了详细的说明，但是，所属技术领域的技术人员能够理解，在不脱离本发明宗旨的前提下，还可以对上述实施例中的各个具体参数进行变更，形成多个具体的实施例，均为本发明的常见变化范围，在此不再一一详述。

Claims (4)

1.一种基于横观各向同性岩体三维蠕变模型的参数辨识方法，包括以下步骤：
（1）将岩石试样制成层理水平岩样和层理垂直岩样，进行横观各向同性岩体的三轴或单轴压缩蠕变试验；
（2）采用陈氏加载法进行处理步骤（1）的试验结果，得层理水平岩样和层理垂直岩样的分级蠕变试验曲线；
（3）三维空间中，设xoy为层理面所在平面，采用以下层状岩体的三维蠕变本构方程，进行相应应变计算：
Figure DEST_PATH_IMAGE002
——式（I）；
Figure DEST_PATH_IMAGE004
——式（II）；
Figure DEST_PATH_IMAGE006
——式（III）；
Figure DEST_PATH_IMAGE008
——式（IV）；
Figure DEST_PATH_IMAGE010
——式（V）；
Figure DEST_PATH_IMAGE012
——式（VI）；
式（I）至式（VI）中：ε x 为x方向的正应变；ε y 为y方向的正应变；ε z 为z方向的正应变；γ yz 为yoz平面的剪应变；γ xz 为xoz平面的剪应变；γ xy 为xoy平面的剪应变；μ hh 、μ hv 分别为平行和垂直层理方向的泊松比；Εh M、ηh M、Εh K、ηh K，及对应的n 1 、n 2 、n 3 、n 4 分别为垂直层理面方向的独立拉压蠕变参数；t为时间；σ x 为方向的正应力，σ y 为方向的正应力，σ z 为方向的正应力；
（4）对于层理水平岩样，xoy平面为各向同性平面，则σ x= σ y= σ 3 ，σ z= σ 1 ，ε z= ε0轴，其应变采用式（III）计算，式（III）变换为式（VII）：
Figure DEST_PATH_IMAGE014
——式（VII）；
对于层理垂直岩样，σ x= σ y= σ 3 ，σ x= σ 1 ，ε x= ε90轴，其应变采用式（I）计算，式（I）变换写为式（VIII）：
Figure DEST_PATH_IMAGE016
——式（VIII）；
式中，σ 1 为轴向大主应力；σ 3 为环向小主应力；
（5）基于式（VII），结合层理水平岩样的蠕变试验曲线，通过Levenberg-Marquardt非线性最小二乘法进行拟合得到参数组合Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 ；
（6）若层理水平岩样和垂直岩样采用相同的应力水平进行蠕变试验，则将所述Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 代入式（VIII）中；根据层理垂直岩样的蠕变试验曲线结果，拟合得到Εh M、ηh M、Εh K和ηh K；
（7）将所述Εh Mn 1 、ηh Mn 2 、Εh Kn 3 、ηh Kn 4 和所述Εh M、ηh M、Εh K、ηh K的计算结果相除，得出n 1 、n 2 、n 3 和n 4 ，得横观各向同性岩体的拉压蠕变参数。
2.依据权利要求1所述的参数辨识方法，其特征在于，在所述步骤（6）中，若层理水平岩样和垂直岩样采用不同的应力水平进行蠕变试验，需根据层理水平岩样的试验结果，通过回归分析计算得到层理垂直岩样应力水平对应的Εh Mn 1 、ηh Mn 2 、Εh Kn 3 和ηh Kn 4 。
3.依据权利要求2所述的参数辨识方法，其特征在于，所述回归分析中回归曲线采用非线性方程。
4.依据权利要求3所述的参数辨识方法，其特征在于，所述非线性方程为指数方程、幂函数方程、对数方程中的任一种。
CN201910338061.1A 2019-04-25 2019-04-25 基于横观各向同性岩体三维蠕变模型的参数辨识方法 Active CN110174307B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910338061.1A CN110174307B (zh) 2019-04-25 2019-04-25 基于横观各向同性岩体三维蠕变模型的参数辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910338061.1A CN110174307B (zh) 2019-04-25 2019-04-25 基于横观各向同性岩体三维蠕变模型的参数辨识方法

Publications (2)

Publication Number Publication Date
CN110174307A CN110174307A (zh) 2019-08-27
CN110174307B true CN110174307B (zh) 2021-06-22

Family

ID=67690116

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910338061.1A Active CN110174307B (zh) 2019-04-25 2019-04-25 基于横观各向同性岩体三维蠕变模型的参数辨识方法

Country Status (1)

Country Link
CN (1) CN110174307B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110793858B (zh) * 2019-11-19 2022-07-15 贵州大学 一种测定岩体中层理力学参数的方法
CN114112651A (zh) * 2020-08-27 2022-03-01 中国石油化工股份有限公司 一种用于人造岩心的岩石动静态力学参数转换方法及系统
CN113836619B (zh) * 2021-09-09 2024-01-30 中国铁路设计集团有限公司 一种拟合任意徐变系数曲线的方法
CN114894616B (zh) * 2022-04-15 2023-06-06 安徽理工大学 一种基于变形模量的岩石蠕变模型粘、弹性参数获取方法
CN114838852B (zh) * 2022-05-10 2023-05-30 重庆科技学院 一种确定地质应力场方向的实验装置及其实验方法
CN115859836B (zh) * 2023-02-22 2023-05-02 西南交通大学 一种岩体蠕变本构模型参数识别方法、装置、设备及介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103995957A (zh) * 2014-04-18 2014-08-20 华东理工大学 针对蠕变损伤及体积型缺陷的承压结构的安全评定方法
CN104020276A (zh) * 2014-06-24 2014-09-03 中国石油大学(北京) 横观各向同性页岩储层岩石力学参数的确定方法
CN105571932A (zh) * 2015-12-17 2016-05-11 河海大学 一种沥青混合料横观各向同性的测试方法
CN106484936A (zh) * 2015-09-02 2017-03-08 中航商用航空发动机有限责任公司 一种考虑应力松弛效应的高温部件的蠕变疲劳寿命的评估方法及装置
CN106769482A (zh) * 2017-03-10 2017-05-31 中国石油天然气股份有限公司 一种静态横观各向同性岩石力学的实验和分析方法
CN106932275A (zh) * 2017-03-10 2017-07-07 中国石油天然气股份有限公司 一种横观各向同性岩石力学测试和分析方法
CN108931448A (zh) * 2018-05-07 2018-12-04 华南理工大学 一种高铬钢材料热力学响应及疲劳-蠕变损伤的预测方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103995957A (zh) * 2014-04-18 2014-08-20 华东理工大学 针对蠕变损伤及体积型缺陷的承压结构的安全评定方法
CN104020276A (zh) * 2014-06-24 2014-09-03 中国石油大学(北京) 横观各向同性页岩储层岩石力学参数的确定方法
CN106484936A (zh) * 2015-09-02 2017-03-08 中航商用航空发动机有限责任公司 一种考虑应力松弛效应的高温部件的蠕变疲劳寿命的评估方法及装置
CN105571932A (zh) * 2015-12-17 2016-05-11 河海大学 一种沥青混合料横观各向同性的测试方法
CN106769482A (zh) * 2017-03-10 2017-05-31 中国石油天然气股份有限公司 一种静态横观各向同性岩石力学的实验和分析方法
CN106932275A (zh) * 2017-03-10 2017-07-07 中国石油天然气股份有限公司 一种横观各向同性岩石力学测试和分析方法
CN108931448A (zh) * 2018-05-07 2018-12-04 华南理工大学 一种高铬钢材料热力学响应及疲劳-蠕变损伤的预测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"横观各向同性岩石弹塑性本构模型与参数求解方法研究";王者超等;《岩土工程学报》;20180831;第40卷(第8期);全文 *
"横观各向同性岩石蠕变性质与本构模型研究";王者超等;《岩土工程学报》;20180731;第40卷(第7期);全文 *
"隧道各向异性软岩力学参数反演及蠕变研究";张海洋等;《地下空间与工程学报》;20180831;第14卷;全文 *

Also Published As

Publication number Publication date
CN110174307A (zh) 2019-08-27

Similar Documents

Publication Publication Date Title
CN110174307B (zh) 基于横观各向同性岩体三维蠕变模型的参数辨识方法
CN110176069B (zh) 横观各向同性岩体三维蠕变特性的表征及其表征模型的构建
Li et al. A numerical investigation of the hydraulic fracturing behaviour of conglomerate in Glutenite formation
Zhang et al. Experimental investigation on the influence of high pressure and high temperature on the mechanical properties of deep reservoir rocks
Shan et al. Study on the triaxial unloading creep mechanical properties and damage constitutive model of red sandstone containing a single ice-filled flaw
Yang et al. Numerical study on failure behavior of brittle rock specimen containing pre-existing combined flaws under different confining pressure
Shirani Faradonbeh et al. Fatigue failure characteristics of sandstone under different confining pressures
Ma et al. Simulating strength parameters and size effect of stochastic jointed rock mass using DEM method
Ma et al. Influence of pore pressure on tensile fracture growth in rocks: a new explanation based on numerical testing
Han et al. A comprehensive investigation on the shear properties of interlayer shear weakness zones
Gui et al. Numerical study of the circular opening effect on mechanical behaviour of rock under confinement
Yang et al. Influence mechanism of coal thickness effect on strength and failure mode of coal–rock combination under uniaxial compression
Fu et al. A damage hypoplasticity constitutive model for cemented sand and gravel materials
Eremin Numerical simulation of failure of sandstone specimens utilizing the finite-difference continuous damage mechanics approach
Eremin Influence of the porosity on the uniaxial compressive strength of sandstone samples
Mohamed et al. Uniaxial compressive strength of composite rock material with respect to shale thickness ratio and moisture content
Li et al. Characterization of Marcellus Shale fracture properties through size effect tests and Computations
He et al. Analytical and numerical solutions for shear mechanical behaviors of structural plane
Baczynski Hoek–Brown rock mass: Adjusting geological strength index for Directional Strength
Mardalizad et al. An investigation in constitutive models for damage simulation of rock material
Wang et al. Effects of differences in mechanical parameters of media on mechanical properties and failure form of composite samples
Zhou et al. Effect of bedding orientation and spatial variability of stratification shear strength on stability of transversely isotropic rock slope
Packulak et al. Practical aspects of boundary condition selection on direct shear laboratory tests
Arzúa et al. Study of the mechanical behavior of intact and jointed rocks in laboratory with particular emphasis on dilatancy
Zhang et al. Stability analysis of channel slope based on FEM strength reduction

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